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ABSTRACT: 

This report contains further considerations of the slip parameter for 
unsteady flow. Several important gas/particle combinations are represented 
and conclusions regarding electrical efficiency are given. Also, a computer 
study which solves the combined Laplace /Poisson equations for geometries 
representing the EHD channel is included. This program is used to predict 
the effect of voltage- scheduling electrodes, and the predictions are checked 
out experimentally. 
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This is the final year-end report under the EHD research contract. The 
work was sponsored by the Naval Air Systems Command under the cognizance of 
Dr. H. R. Rosenwasser. 

One Engineers thesis was generated during this period and it is included 
in the Appendix. A paper was presented at the 7th Intersociety Energy Conver- 
sion Engineering Conference held in San Diego, California. A progress paper 
was given also at the energy conversion program review held at China Lake, 
California. A paper entitled, "Charge Depletion Mechanism Operative in Space 
Charge Flows," coauthored with LCDR Bohley has been accepted for publication 
in the Journal of Applied Physics. 
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I . INTRODUCTION 



EHD devices should have stage efficiencies above 10% if they are to 
achieve a competitive standing as electrical power generators. Although the 
highest efficiencies are calculated for low mobility particles, a low slip 
parameter for steady flow may yield incorrect efficiency values if there 
exist significant levels of turbulence in the flow. This unsteadiness may 
arise from several sources, a particular one being the mixing region of two 
flows. Mixing of jets is designed to improve diffuser performance and hence 
the EHD generator efficiency.^ Another important source of turbulence in a 
channel is the boundary layer; the generator is inherently a high surface-to- 
volume ratio device and it is expected that the boundary layers are fully 
developed. Turbulent boundary layers would then introduce significant levels 
of unsteadiness. 

2 3 

As shown in previous reports 5 , a consideration of unsteady flow has 

led to the definition of a more complete slip parameter. The important equations 

are summarized on page 2. The momentum equation used is a simplified one and 

4 

more exact descriptions would seem desirable; however, Barreto has shown that 
major effects are included in the above mentioned description. Certainly, 

r ^ y 

many treatises are today available^’ ’ which describe the drag of micron-sized 

7 

particles. Boothroyd , for example, in a recently published book summarizes 
some of the information. It is shown in Boothroyd' s book that for Reynolds 
numbers below ten and for particle-to-fluid density ratios greater than 1000, 
the equilibration time given by Stokes law is quite accurate. Another criticism 
of the analysis may be that the perturbation approach is only good for small 
efficiencies. This is indeed true; however, there is no compromise needed 
for the anemometer application described in previous reports. The efficient 
generator would have to be described with the full complement of conservation 
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PARTICLE MOMENTUM EQUATION 
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m = mass of particle 
q = charge of particle 
H = viscosity of fluid 
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equations . One last area of error might be the fact that modeling spherical 
particles is not altogether realistic since droplets may easily deform and 
since information on solids indicates a variety of shapes. Here, however, a 
correction factor may be all that is needed to remedy the situation. 

Perhaps the toughest problem in EHD is the knowledge of the size and 
charge distribution of the charged particles. Some information as to stability 
of droplet size and stability of charge per particle is also needed. For example, 
droplets may break up or coalesce under the pressures that exist under power 

Q 

generation conditions. Some charge depletion may occur in the interelectrode 
space. Hence, only after an injector’s operating characteristics are accurately 
known can improvements be incorporated toward optimum conditions . A tough pro- 
blem also exists in the open-cycle EHD channel if solid particles enter the 
atmosphere in large quantities because this represents a health hazard. 

In this report further consideration is given to the slip parameter for 
unsteady flow. Also, a computer study and some experimentation relevant to 
potential distributions inside the generator are included in the Appendix. 

There, some preliminary work on the measurement of mobility by the Rutherford 

o 

method as outlined by Darrow can be found. The computer study was stimulated 
by the work of Minardi"^. 
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II. SLIP PARAMETER 



The criterion previously stated for slip in turbulent flow is 



V’ (u>) 
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( 1 ) 



Now for particles of density , the mass is given by 



m = 5 <3 

p 



(?" 



r3 P, 



( 2 ) 



where ^ is a correction factor which accounts for non-spherical shapes. 
This factor can best be defined empirically. 

Hence 
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Let $ = the effective particle radius 



~2 

2 p p R 

t v“ = 9^-“ 



(4) 



Figure 1 shows two sets of curves bracketing a region for which the con- 

-3 -2 

tribution of t v iu to the slip is between 10 and 10 . Water and mercury 

droplets are shown because they have been considered most extensively, but 
solid particles may be advantageous for some applications. Polystyrene 
latexes are commercially available in average diameters of 0.1 pm to 100 pm. 
Returning to Figure 1, for less dense particles the curves move to the right 
which is a desirable effect from practical considerations since it is easier 
to obtain and charge larger diameter particles. The curves in Figure 1 can, 
of course, be generalized to any particle/gas combination by plotting 



p co/p, vs. i^R for whatever values of slip are desired. 

r J) ' ~CO 
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Frequency Response (Hz) 




Particle Radius (meters) 



FIGURE 1. Frequency Response as a 

Function of Effective Radius 
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In order to arrive at the total slip we have to integrate S(tu) . If we 
define 

00 

J V(uj) d(o = U^rms) (5) 

o 

then 

00 

s 5 lr + ir I T v “ da> ( 6 ) 

00 00 O 

One would have to know experimentally the spectral density distribution, 

V(iu) . As reported in the literature'*”*', with increasing frequencies V(io) 

-5/3 -7 

is found to be proportional first to tu and later to u) , which corresponds 
to a rather abrupt drop at the high frequencies. For illustrative purposes, a 
simplification will be made here by assuming that V(tu) is a constant value until 
and zero thereafter (i.e., a step function). The resulting integral becomes 

1 06 

— J t V((u) du> = T y ttXjL (7) 

00 o 



and hence the slip becomes 



S 




00 



(8) 



This is reminiscent of our earlier statements to the effect that considera- 
tion of the highest frequency for which V(tu) has a large intensity is all that 
is necessary in the simplest cases. 

We can now define the electrical efficiency or internal efficiency of 
. 12 

energy conversion 

\ ** uE*E = 1 + S ^ 

1 V-E 
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Upon substitution of Equation 8, the above becomes 







( 10 ) 



Clearly, if kE/U is small compared to one but not t tu , then the efficiency 

GO V _L 

will be low, perhaps falling below the 10$ value alluded to earlier in this 
report . 



III. DISCUSSION 

While it is difficult to talk in generalities, some appropriate comments 
can be made here about the mixing of two jets. When a supersonic jet carrying 
the charged particles is injected into a subsonic flow of a low molecular 
weight gas, a highly efficient diffuser results which improves generator per- 
formance by the high pressure operation^. Suppose we consider a jet of mercury 
at 300 m/sec issuing from a nozzle one mm in exit diameter. The size of the 
largest eddies, which initially carry most of the energy, can be no larger 
than the nozzle diameter. In fact, most of the eddies will be generated at 
the edges of the jet as shown in Figure 2, so that taking the nozzle diameter 

gives the lowest frequency to the fluctuations . An eddy characteristic time 

-6 -5 -1 

is of the order of 3 x 10 sec corresponding to a frequency of 3 x 10 sec 

-8 

Now, returning to Figure 1, we must have droplet radii of about 10 m to 
realize a low total slip for this kind of flow. 

This diameter is quite small and droplets may be unstable. Calculations 
of the drag for this type of particle are complicated by the fact that the 
particles do not interact in the fashion of a continuum if conditions are near 

STP. Moreover, the thermal control of the generator system may be quite com- 

-8 

plex in order to get 10 - m particles from a condensing system (powders, on the 
other hand, would be nearly impossible to handle). 
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FIGURE 2. Mixing of Jets Unsteadiness 



What can we surmise will happen if the particles are substantially greater 

-8 

then the low slip value of 10 m? First of all, the diffuser effect will still 
take place but the particles will not be intimately coupled to the flow. Particle 
inertia will influence their trajectories. This could have a strong effect on 
the growth of the particles in the mixing region as well as on the mixing itself. 
Mixing of charges, it will be recalled, is one of the things that affects break- 
down. Jet profiles may differ when the jet contains a large number of high slip 
particles and the positioning of electrodes and of volt age -scheduling surface is 
thereby complicated. In short, performance estimates based on a steady flow 
slip parameter will differ from actuality and they will differ by underestimating 
losses. 

IV. CONCLUSIONS 

In this report, an analysis has been presented which permits the calculation 
of efficiency when experimental information for the velocity spectral distribu- 
tion is available. The effects of turbulence may vary from larger slip and 
frictional losses to an enhancement of the breakdown potential. The field of 
diagnostics in space-charge flows remains underdeveloped and thus our models 
suffer from the lack of validating experimental information (current -voltage data 
from the generator represents the composite action of many processes and is not 
useful in this regard). 
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ABSTRACT 



In an electrogasdynamic (EGD) generator the flow of a neutral medium 
carrying unipolar electrical charges constitutes an electrical current. 

A model of the charge flow in an EGD generator was constructed for use 
in a computer program. The program is designed to solve Poisson's and 
Laplace's equations for both axisymmetric and two-dimensional geometries. 
Schlieren photographs of the charge cloud were used to determine the 
charge cloud profile required by the program. Computer generated pre- 
dictions agreed with three known solutions to Poisson's equation. 

Computer predictions of the effects of space charge flow modification 
were obtained. Space charge flow was modified both by increasing flow 
speed and by manipulating the space charge electric field. Experimentally, 
this was accomplished by increasing flow stagnation pressure and by 
application of 'a separate controllable electric field. Experimental 
results compare favorably with computer predictions. Some measurements 
were also made of the mobility range of the charged particles. 



2 



TABLE OF CONTENTS 



I. INTRODUCTION - 8 

II. THEORETICAL CONSIDERATIONS 10 

III. PROGRAMMING CONSIDERATIONS 16 

IV. DESCRIPTION OF APPARATUS 21 

V. EXPERIMENTAL PROCEDURE 23 

VI. RESULTS AND DISCUSSION 25 

VII. CONCLUSIONS 29 

VIII. RECOMMENDATIONS 30 

APPENDIX A: NORMALIZATION PROCEDURE 31 

APPENDIX B: FINITE DIFFERENCE EQUATION DERIVATION 33 

APPENDIX C: PROGRAM DESCRIPTION 36 

APPENDIX D: MOBILITY MEASUREMENTS 44 

FIGURES 46 

COMPUTER PROGRAM 65 

BIBLIOGRAPHY 83 

INITIAL DISTRIBUTION LIST 84 

FORM DD 1473 85 



3 



LIST OF ILLUSTRATIONS 



1. EGD Generator Schematic 46 

2. Sample EGD Geometry and Boundary Conditions 47 

3. Parallel Electrode Potential Distribution at Breakdown -- 48 

4. Aerosol Injector and Collector Probe 49 

5. Steam Jet Schlieren Photographs 50 

6. Steam Jet Profile 51 

7. Guard Ring Schematic 52 

8. Corona Current and Convection Current Ratio vs. Corona 

Voltage 53 

9. Potential Distribution, P Q = 8 psig, p g0 = .0103 C/m 3 — 54 

10. Potential Distribution, P Q = 14 psig, p g0 = .0084 C/m 3 -- 55 

11. Centerline Distribution of Slip Parameter 56 

12. Potential Distribution, P e0 =0.0 C/m 3 (LAPLACIAN) 57 

13. Electric Field Distribution, P = 14 psig, 

P eo = .0084 C/m3 2 58 

14. Collector Current Increase Due to Guard Ring 59 

15. Potential Distribution, P = 14 psig. Guard Ring 

Voltage = 0 60 

16. Potential Distribution, P„ = 14 psig. Guard Ring 

Voltage = lkV 2 - 61 

17. Potential Distribution, P = 14 psig, Guard Ring 

Voltage = 2kV - - 62 

18. Comparison Between Laplacian £ field Distribution and 

Poisson E Field Distribution Using the Guard Ring 63 

19. Mobility Measuring Schematic 64 



4 



TABLE OF SYMBOLS AND ABBREVIATIONS 



A Normalized potential; Argon; ampere 

A. . Normalized potential at grid point (i,j) 

* 

O 

B Magnetic flux density (Wb/m ) 

C Constant resulting from normalization; Coulomb 

E" Electric field intensity (V/m) 

EGD Electrogasdynamics 

F a __ Electric field due to voltage applied to physical geometry (V/m) 
app 

E^ Breakdown field strength (V/m) 

E $c Electric field due to space charge presence (V/m) 

E Normalized E field from computer solution 

F Fahrenheit 

HL Mesh interval 

h Q Streamline stagnation enthalpy 

I G Collector current (microamperes) 

Ij Corona current (microamperes) 

_ o 

J Current density (A/m ) 

k Mobility (m^/V-s) 

kV Kilovolts 

M Mach number 

m Meters 

mm Millimeters 

N2 Nitrogen 

P Ambient Pressure (psig) 

P Q Stagnation pressure (psig) 

psig Pounds per square inch gauge 



5 



psia Pounds per square inch absolute 

R Normalized jet radius 

r Dielectric gas jet radius (centimeters) 

r Q Injector nozzle radius (centimeters) 

5 Streamline entropy 

SOR Successive over-relaxation 

V Volts 

Vp Total particle velocity (m/s) 

7^ Dielectric flow velocity (m/s) 

Z Normalized axi symmetric axial coordinate 

y Ratio of specific heats 

6 Slip parameter 

e Q Dielectric constant for vacuum and gas (V-m/A-s) 
p Normalized space charge density 

P e Space Charge density (C/m ) 

p go Space charge density at injector nozzle exit plane (C/m ) 
a Maximum error in approximation to jet profile 

<)> Potential function 

oj Over-relaxation factor , vorticity 



6 



I. INTRODUCTION 



Electrostatic phenomena involve many everyday experiences, such as 
walking on a thick rug and then touching a grounded object. Practical 
applications of electrostatic devices range from Xerography to pollution 
control devices or to the van der Graaf generator [1]. Recently, interest 
in electrical power generation has brought attention to the electrogas- 
dynamic (EGD) generator. This direct energy conversion device is ana- 
logous to the van der Graaf generator and operates by convecting charges 
with a gas against an electric field to a point of collection. The 
charge flow constitutes an electrical current. The charge movement is 
accomplished by an exchange of momentum between the flowing dielectric 
medium and the charges. 

The common denominator of electrostatic devices is the existence of 
free charges in a neutral medium or dielectric -- for the present pur- 
poses a fluid. The flow of a medium with unipolar charges is termed 
space charge flow and the understanding of this flow is critical to the 
design of the EGD generator. This thesis consists of a computer model 
of the space charge flow to be used as an EGD generator design tool. 

The computer model has some limitations inherent to the programming 
itself and these are discussed in Section III. In addition, parameters 
such as charge mobility and velocity as well as the inputs of geometry 
and jet configuration are needed. These are obtained with an experimental 
set-up discussed in Section IV. The program is used to predict regions 
of breakdown, given an electrode geometry and a space charge cloud. In 
addition, the program is used to predict the effect of guard electrodes. 
These electrodes are placed for the purpose of 'smoothing' the electric 
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field in the conversion region, and both computer predictions and 
experimental results are given in Section VI for the effect of the 
guard ring. 
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II. THEORETICAL CONSIDERATIONS 



The EGD process requires producing, transporting and collecting the 
charges which comprise the electrical current. One means by which 
charges may be produced is by a corona discharge. Basically this scheme 
consists of two electrodes, which may be a needle and a concentric ring; 
the concentric ring or attractor electrode, when charged to a high posi- 
tive electrical potential with respect to the grounded needle, induces the 
corona discharge. The corona breaks down a neutral gas initially into 
its electrical components: molecular ions and free electrons. The corona 
needle will attract the positive ions, while negative charges will tend 
to migrate to the attractor ring. The EGD objective, instead, is to 
transport a current of negative charges away from this attractor ring. 

An ion's susceptibility to movement in an electric field is termed 
it's 'mobility', (k). Thus the drift velocity of the ion in the corona 
field is determined by the product of the ion mobility and the electric 
field strength to which the ion is exposed. Since the mobilities of 
molecular ions and free electrons are high, they would migrate rapidly 
to the corona electrodes and be lost to EGD power generation. However, 
when a supersaturated vapor is present, negatively charged droplet 
(micron) sized particles are nucleated which, for a given charge, will 
have much lower mobility than an ion [2,5]. Hence, the charged particle 
is more suited to EGD applications [6]. To move the charged particle 
against an electric field, the particle must be well coupled to the 
dielectric flow. When the dielectric gas transports the charge to a 
downstream collector, against an electric field (E), the total charge 
velocity (Vp) may be described as the sum of the dielectric velocity 
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(v ) and the drift velocity (kE). 

oo 

v p =v ro +kE (1) 

In order to obtain the highest degree of coupling, the drift velocity 
of the charged particle must be low compared to the dielectric velocity 
[4]. The charged particles are carried out of the charge production 
section against an electric field in the conversion section, to the charge 
collector. This collector consists of an equipotential metallic surface. 

A path through a load to ground is provided to complete the circuit 
started at the grounded corona electrode. See Figure (1) for a 
schematic of the EGD generator. 

There are several limitations that pertain to the conversion process. 
For an EGD generator to operate as intended, the neutral gas must remain 
non-conducting. This implies that the electric fields to which the 
dielectric is subjected must remain below the dielectric strength, or 
breakdown field strength, ( E^) , of the gas. Fields above the dielectric 
strength will produce regions of high conductivity through ionization 
[3], The dielectric gas or fluid is subjected to the field resulting 
from the applied potential through which the charged particles are carried 
and to internal fields generated by the presence of the charged particles. 
Breakdown from external fields occurs as a result of the electrode geo- 
metry. For a given voltage, the existence of a small radius of curva- 
ture will produce a large electric field. The field will exceed the 
dielectric strength when the radius of curvature is small enough. The 
range of allowable operating voltages is thus a function of the elctrode 
spacing. In turn, the spacing of the electrodes is also a function of 
the fluid dynamics of the dielectric, of possible charge depletion 
mechanisms, of the operating pressures, and of generator size requirements. 
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Total charge per particle is limited by the balance between the 
particle surface tension and charge repulsion forces, and is known as 
the Rayleigh limit. In practice, however, the total charge per particle 
is governed by the charge production method. For particles charged in 
a corona field, and which grow to micron size, the number of charges is 
approximately two orders of magnitude less than that predicted by the 
Rayleigh limit [6]. These limits must be considered in the modeling of 
the space charge flow in an EGD device. 

In general, the equations needed to describe the EGD flow will include 
the momentum, continuity and energy equations. These would be necessary 
to determine the dielectric jet profile and hence the space charge 
density distribution. However, if the profile can be determined 
experimentally, and if the energy drained from the flow is low, the 
required equations are reduced to three: 

Faraday's Law 



_ -f _ 3 B 

v * E - n 



Gauss' Law 



v 




Charge Conservation 



E_ = electric field (V/m) ~ 

B = magnetic flux density vector (Wb/nr) 



o 

= space charge density (C/m ) 

= dielectric constant for vacuum and gas 
(V-m/A-s) 



v • p v = 0 (2) 

e p 

where the total velocity of the charged particle is given by: 

*p - ». + k < r app + T sc> (3) 

and is a function of the dielectric flow speed, particle mobility, 

applied electric field (E___) , and the space-charge-produced electric 

app 

field (E sc ). 
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For the EGD generator, B = 0, and Faraday's Law takes the form, 
v x E = 0. This allows E to be defined as the negative gradient of some 
potential function <j>, so that E = - v<}>. Substituting into Gauss' Law, 
the resulting equation is known as Poisson's equation, 

? p p 

O = ~ (4) 

e o 

This equation is valid in the two regions of interest within the EGD 
channel. One, the neutral gas jet, contains the charged particles, while 
the second, outside the gas jet, contains no space charged. In the 
space-charge-free region, p g = 0, and Eq. (4) reduces to the more famil- 
iar Laplace equation. Figure (1) depicts the two regions. 

Although the kinetic energy of the injected gas jet is generally 
augmented by the transfer of momentum from a primary flow, this method 
of flow augmentation is not to be used here because of the resulting 
unsteady motion of the jet [7]. Rather, the injected flow is augmented 
as described in Section IV. 

Reference [8] states that the jet profile may be determined by con- 
sidering the equation of conservation of charges. The radius of the gas 
jet (r) may be described as a function of the initial radius (r Q ), 
initial charge density ( p eo ) > the applied electric field and the axial 
coordinate (z). 



This is based on the assumption that the charge density is a function only 
of axial distance from the point of injection and that the particle 
velocity due to coupling with the neutral gas and applied field is much 




(5) 
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greater than that due to the space chrage field, i.e.. 



v + k E » k E 
°° app sc 



( 6 ) 



Neglecting space charge effects disallows any radial velocity. Yet 
there must be some component of velocity. For flow along a streamline, 
Crocco's theorem for a steady state, which relates enthalpy (h Q ) and 
entropy (S) gradients to vorticity, is given by: 



v x v x v =vh - TvS 

00 00 0 

where v x 7 = « is the vorticity. For the streamline defining the 
boundary of the gas jet, such gradients will be very large, giving rise 
to vorticity. The radial velocity component introduced will give a 
larger jet profile than that predicted by Eq. (5), leading to a reduced 
charge density at any given axial distance. 

In general. Equations (2) and (4) are coupled through the p g term, 
and must be solved simultaneously. Expanding the charge conservation 
equation: 



!i|_7 + !£8_ p + p |_7 + 7 !-p =0 

r 3r P r 3r ^e M e 3z p p 3z y e 



For particles of low mobility, and fields below breakdown, v^ will 
predominate, and Eq. (6) will hold. If the jet is assumed to be flowing 
normal to two plane parallel electrodes: 



3 

3r 



i-v 
3r oo 




app 



= 0 



Then an assumption of an initial homogeneous distribution of charge at 
the jet entry plane will remain true at any position downstream. Thus, 



3 

3r p e 



= 0 
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and the charge conservation equation becomes: 

v ’ ( p e v p^ p e 37 v p + v p d~z p e " ^ 

This allows the cloud shape and hence the charge density to be determined 
solely by the fluid dynamics of the dielectric jet, and uncouples the 
charge conservation equation from Poisson's equation. Another approach 
to this problem is to determine the jet profile experimentally. An 
approximate analytical description is obtained by fitting a parabola to 
this profile. This approach is convenient in this work and the 
approximation becomes an integral part of the computer program. 

We can now discuss the assumption of Eq. (6) that the charged par- 
ticles are affected only by the gas jet and the applied field. When 
charged particles are injected into the conversion region, they will face 
the field produced by the charged particles already in the conversion 
section. Eventually, these charged particles will face fields which may 
repel all charges but those with the highest degree of viscous coupling 
to the gas flow. In this manner, space charge effects provide a velocity 
"filter" through which the gas jet must move [9]. 

The expression for current density (J) in the dielectric jet is 
given by: 

o = p e Vp 

This expression would require a greater charge density as space charge 
effects slow down the charged particles in order to maintain a given 
current density. But as mentioned previously on Page 11 » breakdown of the 
dielectric will occur if the space charge electric field strength is 
greater than the dielectric breakdown strength. This tradeoff limits 
the performance of the EGD generator. 
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III. PROGRAMMING CONSIDERATIONS 



The computer program was written for the numerical solution of Eq. 

(4), since closed form solutions to Poisson's equation are difficult to 
obtain, except in very simple cases. The program is capable of solving 
Poisson's equation for either an axisymmetric or two-dimensional geometry. 

To facilitate the solution, Poisson's equation was normalized as 
described in Appendix A. In the axisymmetric form, the equation to be 
solved becomes: 

R A R + A RR + A ZZ = ' Cp ^ 

where A is the normalized potential and R is the normalized radius. The 
term on the right side is a result of the normalization described in 
Appendix A. In the two-dimensional equation, there is a change of 
variables and the radius term drops out to give the form: 

A XX + A YY = ‘ Cp 

By covering the geometry to be studied with a uniform square mesh, 
the problem may be defined in terms of a finite number of discrete points, 
at each of which the solution to Eq. (7) is desired. The boundaries of 
the mesh are determined either by the presence of electrodes at a con- 
stant potential, or by a zero normal derivative of the potential <f>. 

The zero normal derivative is found on an axis of symmetry, on a boundary 
with no normal current flow, such as an insulator, or on a boundary 
between two electrodes in the absence of space charge and end effects. 

End effects result in a distortion of the distributions about the end of 
a pseudo-infinite plate because of the radius of curvature of the end. 
Figure (2) depicts a sample geometry with appropriate boundary conditions. 
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The derivation of the finite difference equations for the two forms 
of Poisson's equation is given in Appendix B. After solving for the 
potential at a particular point, the form of the equation is: 



A i ,J 



hAp ♦ A j+1J + A iJ+1 , 4|<A jjjtl -A^J.,) 



where the subscripts represent mesh row and column, respectively. For 
points governed by a boundary condition of a zero normal derivative, as 
on an axis of symmetry, the general equation becomes: 




= 1.333 A 



i J+1 



.333 A 



i »j+2 



The derivation of this equation is also given in Appendix B. 

The solution method chosen is an iterative process known as "success 
ive over-relaxation" (SOR), or the "extrapolated Liebmann method [10]." 
Over-relaxation introduces a factor (co), by which the solution at the 
point in question, as a result of the previous iteration, is taken into 
account. By assigning the proper value to this factor, convergence rate 
of the solution may be maximized; however, w is limited by the following 



1 < a) < 2 



The method of selecting a value for u> may be found in textbooks covering 
iterative solution methods such as Reference [10]. For this thesis u> is 
a constant equal to 1.5. The form of the resulting equation then 



becomes : 










A n+1 
A i J 


= A n + “ 

A i J 4 


HL 2 Cp + Aj +1 




+ A n+1 
A i , j+1 




+ dk. / A 11 _ ) 

+ 2R \ A i , j+1 A i ,j-l / 


- 4 a ?j 
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This shows that the solution at a point, on the (n+l)st iteration is 
dependent on the solutions at surrounding points obtained on either the 
(n)th or (n+l)st iteration. Convergence is determined by comparing a 
specified value, DELU, to the maximum potential change per iteration. 

The value of DELU is a data item required by the computer program and 
remains constant until changed by the user. 

This simplification of the governing equations introduced in 
Section II allows the space charge density to be described as a function 
of the initial space charge, of jet radius and of the initial jet radius. 
The expression, as shown by Minardi, is: 




This leads to a charge density varying with the longitudinal coordinate 
only and a constant distribution in the plane normal to the dielectric 
flow. The use of a known nozzle radius, i.e., intial jet radius, of a 
gas flow speed, and of a corona current (Iy), allows determination of 
the injected space charge density. 

Also to be considered when modeling the charge distribution is the 
collection process. Charged particles cannot continue to be collected 
indefinitely. Eventually, at some point downstream from the initial 
point of collection, all but a negligible part of the charged particles 
that are to be collected, will have been collected. Further charge 
contributions to the collector current will have a negligible effect on 
the potential distribution. It is necessary to determine where this 
takes place, as the boundary conditions on the problem include a zero 
charge density. Hence, the space charge contained in the dielectric jet 
is not to be carried up to the boundary of the region of interest. 

Using various collection lengths in the computer program, it was found 
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that allowing space charge to exist over 80% of the collector length 
gave results which changed little, less than 3%, as this collection 
length was increased. 

In general, the jet profile does not coincide with a mesh point at 
any given axial coordinate. However, the program is not capable of 
handling a variable mesh size. In order to differentiate clearly 
between the regions of space charge and no space charge, the program 
approximates the gas jet with a series of steps in order to meet the 
constant mesh interval limitation. This is explained with the help of 
the figure below: 




□I 

If Ar > 2 ~ , the gas jet is considered to be passing through the point 
(i,j+l), and if Ar < 75— , through the point (i,j). The charge density is 
computed according to Eq. (8), giving a value of charge density which is 
then assigned to each mesh point out to, and including, the point on the 
approximated jet boundary. This introduces some error into the calcula- 
tion process. The approximation scheme limits the displacement of the 
dielectric jet boundary to a maximum of HL/2. Then the maximum error 

<°max> is 9' ven bj,: 



a 

max 



HL/2R 
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This would imply that the maximum error would occur at the minimum jet 
radius, the injector nozzle. By selecting the appropriate mesh interval, 
the radius at the nozzle can be made to coincide with a mesh point. The 
minimum radius at which the maximum displacement will occur is then 
downstream from the nozzle. o max = 5% is typical of the various mesh 
sizes used. 

The validity and accuracy of the program was determined by compari- 
son with the known solution for several problems. A geometry consisting 
of finite parallel plane electrodes, approximating infinite parallel 
plane electrodes, with a gap spacing of 2 cm and a potential difference 
of 15 kV was used to determine if the program would predict the exist- 
ence of electrical breakdown across the gap as specified in [11]. 

Because of the method of normalization used, as described in Appendix A, 
the result of entering this configuration into the program is a potential 
distribution which varies linearly from zero to one, and a field strength 
of one throughout the gap, indicating that breakdown of the air gap is 
predicted. To attain this solution, it was found that electrode end 
effects had to be minimized, or a very distorted distribution would 
result. A pseudo-infinite plane electrode minimum length to air gap 
ratio of one was required to minimize the distortion. Figure (3) depicts 
the correct potential distribution calculated by the program once the 
paral lei plane electrodes were properly modeled. 

Results of the program were also compared with two results previously 
published by Minardi [2,8], concerning space charge flow between parallel 
infinite plane electrodes; values of the input data used were the same. 

A comparison of the graphical output appeared to be equivalent for the 
two geometries compared. 
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IV. DESCRIPTION OF APPARATUS 



The need to input certain parameters to the program dictated the type 
and extent of the experimental work to be accomplished. This work also 
had to be relevant to the EGD facility that was under operation at NPS 
[4,5]. The program requires both the electrode geometry of interest and 
the dielectric gas jet profile, as inputs, so that charge density distri- 
bution may be calculated. The gas jet profile is also the means of 
modeling the change in a variable such as the dielectric flow speed. The 
effect of an increase in 7 ro was of particular interest because a larger 
7^ would increase the dominance of the first term on the right of Eq. (1). 
Also of interest was the effect of the guard ring electrode on the space 
charge flow. These two experiments, increasing 7 ot and the use of the 
guard ring, were also required for correlation of the trends in experimental 
results with those trends predicted by the computer program. 

The EGD channel to be modeled consists of a charged aerosol injector 
and a collector probe. See Figure (4) for dimensions. The vapor for the 
aerosol is generated by a modified Hotshot model MB-31 electric steam 
boiler, pressure fed with distilled water. To obtain best results, it 
was found necessary to ensure that the entire steam generating system 
was free of any impurities. The corona needle is grounded through a 
Simpson microammeter, while the attractor ring is powered by a Sorenson 
high voltage supply, with a range of zero to ten kV. The power supply 
output is monitored by a Sensitive Research high impedance voltmeter 
connected across the corona attractor ring gap. The collector needle 
is grounded through a Simpson microammeter in series with a diode to 
prevent reverse current flow [7]. 
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In order to determine the jet profile, Schlieren photographs of the 
steam jet were taken. The setup consists of mirrors, a collimated light 
source, a knife-edge and a camera. Photographs of the jet were taken 
with stagnation pressures of 8, 10 and 14 psig. A horizontal knife-edge 
was used to more clearly define the vertical gradient across the jet 
boundary. The Schlieren photographs of the 8 and 14 psig jet are shown 
in Figure (5). The photographs show a high pressure jet (14 psig) 
approximately 100% larger than the given by Eq. (5). The resulting jet 
profiles are shown in Figure (6). 

In order to determine the effects of manipulating the electric field 
distribution, a field entirely separate from that of the corona was 
applied by means of a 'guard electrode', in the form of a ring of stain- 
less steel. Internal diameters used ranged from 3/8 of an inch to 3/16 
of an inch. The ring was mounted on a plexiglass insulator, supported 
by a stainless steel rod. The rod was mounted on a vernier allowing the 
position of the ring to be varied along the longitudinal axis of the 
collector needle. The voltage applied to the guard ring was supplied 
by a Spellman high voltage power supply. The output voltage of the power 
supply was determined by calibration of the power supply voltmeter with 
the Sensitive Research voltmeter. This was necessary because the 
Spellman voltmeter does not indicate line voltage. The guard ring 
schematic is shown in Figure (7). 
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V. EXPERIMENTAL PROCEDURE 



Initial experiments with the guard ring were carried out at eight 
psig. The temperatures used corresponded to vapor states at either 
side of saturated vapor, similar to previous EGD work [4,5,7]. At this 
low pressure, a severe deterioration of the 'convection current' ratio, 
(collector current/corona current), was encountered as ring voltage was 
increased. This appeared to be the result of too low a dielectric 
velocity and hence a high slip parameter (i.e., a high ratio of the drift 
velocity to the dielectric velocity). The charged particles were diverted 
from the col lector, being attracted to the guard ring. For this reason 
and for the desire to improve the generator performance without the ring, 
it was decided to increase the flow rate of the dielectric by increasing 
the reservoir pressure. Stagnation conditions were changed by the 
introduction of a separate 'dry' gas into the steam boiler. 

The baseline characteristics of the EGD channel were determined by 
supplying steam to the corona discharge at eight psig and 240-245°F. 

Corona voltage was varied over a range of zero to three kV, and values 
of corona current and collector current were recorded. The collector 
probe was positioned 5 mm downstream from the nozzle exit plane. Dry 
gases used were both N£ and A. The choice of gas appeared to make no 
difference in the results. Introducing the dry gases was done by 
regulating the flow of gas from a standard high pressure gas bottle. 

The gas was introduced into the boiler through acheck valve operating 
at 10 psig. With the boiler operating at 8 psig, the high pressure gas 
was used to raise the operating pressure of the boiler to approximately 
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13.5 psig. After allowing the rate of steam production to stabilize, 
resultant stagnation pressure was 14 psig. 

To determine the effects of the 'guard electrode', a baseline 
performance with no guard ring was run at 14 psig. Once this was 
established, the trailing edge of the guard ring was placed in one of 
three positions, either 2 mm upstream, coincident with, or 2 mm down- 
stream of the collector needle tip. At each guard ring position, the 
applied potential was varied over a range of 0 - 3.5 kV, and effects on 
collector current were recorded. Rings of various sizes were used. 

When rings smaller than 3/8 of an inch diameter were used, breakdown 
occurred, apparently caused by stray droplets providing a direct path 
for current transmission. Condensation of the steam inside the ring 
was also observed, and was a possible factor contributing to breakdown. 
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VI. RESULTS AND DISCUSSION 



The introduction of a high pressure gas into the system resulted in 
collected currents of two to three times greater than the low pressure 
steam. Figure (8) shows the baseline (8 psig) corona current-voltage 
curve and the corona current-voltage curve for the high pressure steam 
jet. The two curves are representative of the range of corona operation. 
Figure (8) also shows the variation of the value of the 'convection 
current' ratio. 

In the expression: 



+ k E ,„n + k E or 

app sc 



the term v ro may be given by the following expression from Shapiro [12]. 




where : 




oo 



For the pressures considered, 8 psig and 14 psig, the respective Mach 
numbers are 0.832 and 1.02. These correspond to flow speeds, at the 
injector nozzle exit plane, of 402 m/s at P Q = 22.7 psia, and 483 m/s 
at P Q = 28.7 psia. From the Schlieren photographs the profile of the 
jet was determined. Assuming continuity, this leads to a prediction of 
the variation of flow speed in the flow direction. 

The increased effectiveness of the high pressure jet may be shown 
by recalling the slip parameter, 6, defined to be the ratio of the par- 
ticle drift velocity, kr, to the steam jet flow speed. In terms of the 
non-dimensional electric field output of the computer program, E, 



25 



the slip parameter is: 



KE b |E| 




For a jet with particles well coupled to the flow, 6 will be less than 
one. To determine 6, the value of E may be solved for by the computer 
program. From Eq. (8), the value of p eQ , the initial charge density, may 
be determined. For the low and high pressure jets, this results in 

3 

0.0102 and 0.0084 C/m , respectively. Using this density as the program 
input parameter RHOZRO, an |E| field distribution may be obtained. This 
was done for both low and high pressure cases. Other significant para- 
meters used were attractor voltage of 2300 volts, and breakdown field 
strength of 3 x 10^ V/m. Figures (9) and (10) are plots of the equi- 
potential surfaces resulting from the computer solution of the low and 
high pressure jets. Using these results, the variation of 5 along the 

centerline of the EGD channel was obtained, and is plotted in Figure (11). 

-5 2 

The values of 5 are based on an assumed mobility of 4 x 10 m /V-s, 
a value representative of the range of mobilities of the EGD charged 
particle [7]. 

Figure (11) shows that everywhere along the centerline of the EGD 
channel, with the exception of a single point, the high pressure jet has 
a lower slip than the low pressure jet. The single point that differs 
may well be the result of the inaccuracies of the numerical solution. 

The use of a higher value of mobility would tend to give values of 6 
greater than one. However, if the conjectured distribution of [5] is 
approximately correct, particles of higher mobility are a small part of 
the total number of charged particles. 
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From Eq. (3), the total E field has two contributions. One, E , 

app 

is due entirely to the geometry and voltages applied to various parts 
of the E6D geometry. This can be changed only by a physical change in 
the geometry or applied voltages. The potential and £ field distribution 
resulting from the geometry alone may be solved for by Laplace's equation. 
The Laplacian solution represents the potential or F field distribution 
unaffected by the presence of space charge. Figure (10) is the solution 
to Poisson's equation, and it shows the added distortion of the dis- 
tribution due to the presence of space charge, when compared to the 
Laplacian solution. Figure (12). The added distortion appears in Eq. 

(3) as the second contribution, F sc , the space charge field. Figure (13) 
depicts the |E| field distribution corresponding to the potential 
distribution of Figure (10). Since the space charge is of the same 
polarity as the charged particles, their movement is opposed. Elimina- 
tion or reduction of F sc would enhance charged particle velocity, 
making the collection easier. Application of a potential to the guard 
ring is designed to impose an additional external field on the flow. 
However, if the field overcompensates for the space charge field, and 
the ring becomes a dominant sink for the charged particles, the purpose 
of the ring is defeated. 

Figure (14) shows the increase in collector current as a result of 

l 

using the guard ring. The results are shown over a range of voltages 
at the three positions used. It appears that the position achieving 
the greatest increase is that of the trailing edge plane coincident with 
the tip of the collector needle. From Figure (10), this is also the 
area of greatest space charge distortion, implying that the effects of 
the ring are concentrated in a small area. 
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Figures (15), (16) and (17) are graphical computer outputs of 
potential distribution for a grounded ring, and potentials of 1 and 2 kV. 
The experimental results also show the proper trend in collector current 
as ring voltage is increased. That is, the current slowly increases to 
a broad maximum, and then drops off with further voltage increases, the 
maximum occuring at about 1 kV. 

The computer predicted results of using the guard ring is shown in 
Figure (18) as a plot of centerline Laplacian E minus the Poisson E 
field. As this difference approaches zero, the particular guard ring 
distribution approaches the optimum Laplacian distribution. The computer 
results indicate that the E field is definitely modified by the presence 
of the guard ring. The predicted optimum ring voltage, corresponding to 
maximum collector current, appears to lie between ground and 1 kV. The 
point common to all curves lies at the midpoint of the ring, and is prob- 
ably due to the uniformity of the guard ring electric field at this 
point. 
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VII. CONCLUSIONS 



The computer program appears to present an accurate picture of the 
potential and electric field distributions within an EGD device. This 
is borne out by the computer prediction of experimentally observed trends 
in EGD channel performance. The 'convection current' ratio, Ig/I-p was 
found to be affected by the steam reservoir conditions. In addition, 
computer and experimental results agree and show that increased EGD 
performance can be attained by suitable manipulation of the electric 
field distribution, and hence of the space charge flow. 

The agreement found between analytical and experimental results 
suggests that the use of the computer program as a design tool, par- 
ticularly when modeling the space charge flow in a high pressure jet, 
is feasible. 

To be completely general in approach, some further account should 
be taken of particle mobility. Since the drift velocity appears to be 
negligible at high pressures, this is especially true when working with 
particles of high mobility, and/or dielectric jets at low speeds. In 
this case, the slip parameter could go well above 1, and the present 
assumption of zero slip will give completely invalid answers. Some 
preliminary work on determining the mobility range of the charged 
particles used in these experiments is shown in Appendix D. 
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VIII. RECOMMENDATIONS 



In addition to the measurement of mobility, the difference equations 
should be modified to allow for a variable mesh size. Equations of this 
type may be found in a textbook such as Reference [10]. This will allow 
a greater generality in selecting geometries to be studied, and will also 
do away with a great deal of the error involved in discrete jet profile 
approximation. It will almost certainly mean greater computer storage 
requi rements . 

As a design tool, the program can be used to determine the optimum 
guard ring geometry and voltage. From the curves of Figure (18), the 
solutions for the various voltages diverge from the optimum in the 
vicinity of the collector needle. A guard ring geometry capable of 
affecting this divergence, yet maintaining the distribution of the 
remainder of the 1 kV curve, would appear to be an improvement. 

It should also be possible to use the program as a design tool to 
optimize both injector nozzle and collector geometries. A nozzle 
geometry allowing for the use of smaller guard ring sizes may allow 
for more effective guard ring performance, yet forestall breakdown. 

An optimum collector geometry may exist which will allow for greater 
charge collection efficiency resulting simply from the geometry, rather 
than by flow manipulation to achieve the higher collector currents. 

Other applications outside the EGD field may also exist. For 
instance, the program may be used to study the region of space charge 
known as an electrode 'sheath' which divides a neutral plasma and an 
electrode. Also, Reference [4] covers a study of EGD control of separated 
flow. This study may be extended through use of the program to determine 

optimum electrode geometries for control . 
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APPENDIX A: NORMALIZATION PROCEDURE 



The vector form for Poisson's equation is: 



v 
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<f> = - 
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e 



e 
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After expanding in cylindrical coordinates: 



By defining: 
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where: 



h = characteristic channel length of EGD geometry 
V = maximum potential 
and substituting into the various terms: 
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By combining and clearing the left side: 
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The normalized equation becomes: 



R A R + A RR + A ZZ 
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Similarly for the two-dimensional form of Poisson's equation, 
defining: 



X = 



x 

h 




By 



The two dimensional form becomes: 



A XX + A YY = ' CP 

The electric field may be normalized by defining: 

where: E e the normalized value of E 

E^ = the dielectric breakdown field strength 

from: E = -v<j> 



then : 



hE 
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APPENDIX B: FINITE DIFFERENCE EQUATION DERIVATION 



J 

A 



^ I 



i-1 i j+1 



i .j 




— HL 
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i+1 

j+1 



i+1 

j 



i+1 

j-1 



Knowing the value at any point in the mesh, the value at another 

point, nearby, may be found by means of a Taylor series. By writing the 

series in the J direction, expressions for A. and A. • ■, may be 

1 J J I *1 9 J " I 

found. These are: 



= A ij + HLfl ij + 2T iA i: j + 3T iA i;i + 



A. • -i = A, . - HLA! . + tjt— A ! 1 . - A ! 1 + . 



0) 

( 2 ) 



Solving for A! ., neglecting higher order terms where primes indicate 
the derivative with respect to the direction of interest: 



A 1 

A i , J 



- A i . J+1 " A f .J-1 
2HL 



(3) 



To determine the finite difference equation for the second partial 
derivative, add (1) and (2). This results in: 



A 



1 .j+1 + A i ,j-l 



1 >J 
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Solving for A!'-, again neglecting higher order terms: 
I 



A 1 ' 

A i J 



A i,j+1 + A i,j-1 ' 2 A i,j 

ii? 



( 4 ) 



By writing the series in the I direction, solving for values of A^-j ^ 

and A- . , and adding the resulting equations, the second derivative 
1 " I 5 J 

in the I direction may be determined. The resulting equation is: 




A i +1 ,j + A i-1,j 



HL 




(5) 



By combining (3), (4), and (5), and solving for A. ., a finite difference 

1 sJ 

form of Poisson's equation in axi symmetric form is obtained. The 
axi symmetric form is: 



A i , j = 4 



HL 



' Cp + A i+l,j + A i-1 ,j + A i,j+1 + A i,j-1 



+ 2R (A i ,j+l ' A i J-1* 



( 6 ) 



For points along a boundary where the normal derivative is equal to 
zero, the equations for the value on the boundary may be derived in the 
same manner and cast in the form of second order forward or backward 
difference equations. Assuming (i,j-l) is on the boundary of interest, 
the series equations become: 



A. = A, . , + 2 HLA! , , + 
i ,J+1 i,J-l i.J-1 
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(8) 



Multiplying (8) by 4 and subtracting from (7): 
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Neglecting higher order terms: 






= A i " 4 A i , J + 3 A i , J-l 
- 2 HL 



(9) 



Setting the normal derivative equal to zero, and solving for A. • 

1 5 J “ * 



A i , j-l 1 ,333 A i ,j " * 333 A i ,j+l 



( 10 ) 



The successive over-relaxation method (SOR) uses previously calculated 
values to speed up convergence. For a point inside the mesh, the final 
difference equation becomes: 
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For a point on a boundary of zero normal derivative, the equation 
becomes : 
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APPENDIX C: PROGRAM DESCRIPTION 



The computer program is designed to solve a finite difference form 
of Poisson's equation, either in axisymmetric or two-dimensional form, 
at each node defined by a uniform square mesh covering the geometry to 
be studied. In order to do this, the program must be able to do the 
following: 

1. Determine which finite difference equation is to be solved. 

This is a function of the geometry and the node position within 
the mesh. 

2. Determine and store charge distribution at each node. 

3. Solve the difference equation and store the answer. 

4. Determine when a solution to the problem has been attained. 

5. Provide the results in an easily interpreted form. 

The program is made up of two major parts: preparation and calcula- 
tion, and output of results. In general at each point in the EGD channel 
defined by the uniform mesh, there exists a charge density and a potential. 
The charge density may be zero. This can be the case inside of the jet, 
as well as outside. The voltage may either be the result of an electrode 
presence or a potential surface between electrodes. If the voltage is 
due to an electrode presence, this value will remain constant, otherwise 
it will vary as a function of the charge density. The program must be 
able to make this differentiation. In all, the program must know the 
potential at a point, the charge density, and whether the point defines 
a portion of an electrode. In order to simplify the program coding, the 
program stores all values of potential in one array, "A". The values of 
charge density are stored in a separate array, "G" , of the same size. 
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This array is also used to store all information required to determine 
whether the potential is that of an electrode, the potential description 
value (PDV), or is to be changed as a function of the charge density. 

This can be done since a point of constant potential will be an electrode, 
and a charge density cannot exist at the same point. 

The program determines convergence of the solution method by compar- 
ing answers from the (n+l)st iteration with those obtained on the (n)th 
iteration. This will reveal the maximum change in potential as a result 
of the iteration. If the maximum is less than some designated value, the 
iterative procedure is stopped. This method requires storage for the 
answers generated by the (n)th iteration. These are stored in array "B". 
The value of the potential resulting from the previous solution is required 
for the SOR method. This value is supplied from array "B". Hence the 
variables required for the solution of the finite difference equation 
are taken from three separate sources. A superposition of arrays "A" 
and "B" describes the complete situation at any node within the EGD 
channel . 

The program has been written for use on an IBM 360/67 Computer 
System. The output is in both printed and plotted form. The plotting 
was accomplished on the CALCOMP Model 765. The program consists of 
the following subroutines: 

SUBROUTINE NAME PURPOSE 



MAIN (Main Program) 



GSPRD 

GENFIL 
CAL Cl 



General coordination of preparation 
and calculation phase and output 
phase. 

Calculates cloud radius. Determines 
charge density at each point, and 
stores in "G" array. 

Reads geometry input. 

Calculates potential. Determines E 
field components and modulus. 
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WRITER 



Prints values of potential, E field 
and E field components, determined 
in CALC1 . 

Determines lines of constant potential 
and constant E field to be plotted. 
Orients arrays to coincide with output 
labeling convention. 

General coordination of output phase. 
Coordination of labeling and plotting 
of graph, and preparation for plotting 
results. Modified NPGS Computer 
Library subroutine, which allows the 
superposition of three sets of contours 
on one plotted output. 

All data, except the electrode geometry is read in under NAMELIST 
control. This requires a special format for all cards, but provides for 
a more flexible means of handling the data in that there are no fixed 
fields in which the variables must appear, nor must the arrangement 
coincide with that of the NAMELIST statement. 

The first data card required by MAIN is the card denoting the appear- 
ance of the data group. It has the format, starting in column 2; &PARAM. 
PARAM is the NAMELIST name. The first column of all data cards is ignored. 
Each data item may be separated from the following by a comma. However, 
no comma must appear between the NAMELIST name and first data item. 

More than one data item may appear on a card. Each data item is in one 
of the two forms below. 

1. "Variable name" = FORTRAN constant 

2. "Array name" = a set of FORTRAN constants, separated by constants 
[13]. 

The variables read in are listed below. The first five are explained 
with the help of the representative geometry in Figure ( 2). 



CLVL 

FLOP 

OUTPUT 

CONTUR 
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NAME 



DESCRIPTION 



TYPE 

NATRAD 1*4 

NBLOCK 1*4 



NDLPNT 


1*4 


NDLRAD 


1*4 


NOZPOS 


1*4 


BRKDWN 


R*4 


CONVER 


R*4 


DELU 


R*4 



GAPPOT 


R*4 


IH 


1*4 


IW 


1*4 


NCOLS 


1*4 


NDIMEN = 0 


1*4 


= 1 


1*4 


NGLOSS = 0 


1*4 



Mesh column number of attractor inner 
radius. This is the point defined as 
(6,3). Then the value assigned to 
NATRAD is 3. 

The number of discrete points, lines, 
or blocks of points required to 
completely describe the charge dis- 
tribution or geometry of the problem. 

For the figure shown, the attractor 
ring may be defined by a block of 
points, a line and a single point. 

The corona needle, and collector 
needle may each be defined by two 
lines of points. In addition, if it 
is assumed that the collector needle 
collects all the charge, a block of 
points can be used to define the area 
within the dielectric jet where in no 
charge exists. Thus, for the geometry 
and charge distribution of the problem 
shown, NBLOCK = 8. 

Mesh row position of corona needle 
tip. For the geometry shown, NDLPNT = 4. 
Mesh column number of corona needle 
centerline. For geometry shown, 

NDLRAD = 1. 

Mesh row number of attractor exit 
plane. For the geometry shown, 

NOZPOS = 6. 

Breakdown field strength of dielectric 
gas. (volts/meter) 

Characteristic channel length of EGD 
geometry, (centimeters) 

Value compared with maximum difference 
between iterations. If the maximum 
difference is less than desired value 
of DELL), convergence of the iterative 
scheme is signaled. 

Maximum fixed potential covered by the 
mesh. In general, the attractor will 
be at the maximum potential. 

Desired height of output plot measured 
along longitudinal coordinate (IH s 15) 
Desired width of output plot measured 
along radial coordinate ( IW s 9) 

Number of mesh columns (NCOLS * 61). 
Indicates problem is a two-dimensional 
geometry in constant area channel. 
Indicates problem is axi symmetric. 

No charge density decrease accounted 
for in the conversion section other 
than by charge cloud spreading. 
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NGLOSS = 1 


1*4 


Some loss mechanism (recombination, 
etc.) and cloud spreading will be 
accounted for in determining charge 
cloud densities. 


NL 


1*4 


Number of potential contours desired 
in pi otted output (NL s 20). 


NLI 


1*4 


Number of E field contours desired in 
plotted output (NLI < 20). 


NOLINE = -1 


1*4 


Only potential and E field contours 
will be plotted. 


= 0 


1*4 


Only physical geometry (fixed potential) 
points will be plotted. 


= 1 


1*4 


Physical geometry will be superimposed 
on contour plot. 


NROWS 


1*4 


Number of mesh rows (NROWS s 81) 


RH0ZR0 


R*4 


Charge density at exit plane of 
attractor (coulombs/cubic meter) 



The last card in the NAMELIST signals the end of the list. The format 
is &END. Again the symbol ampersand must appear in column 2 of the data 
card. 

The jet radius and charge density distribution are calculated by 
subroutine GSPRD according to the charge profile given by statement 
numbered 35. In the case of a two-dimensional geometry, only the charge 
distribution is calculated, as the jet is assumed to fill the entire 
two-dimensional channel. The charge density will remain constant through- 
out the flow in the channel, unless some charge loss mechanism is 
introduced, signaled by the value of NGLOSS. The jet profile and charge 
distribution for a geometry with an off center corona needle is also 
determined by GSPRD. 

Subroutine GENFIL is used to input the specific values of potential 
and charge density or PDV. This is accomplished with the help of the 
variable NBLOCK, which is used as a DO-LOOP parameter. The loop causes 
the program to read the exact number of data cards defined by NBLOCK. 

Since each block defined by NBLOCK contains points of common value, the 
entire block can be represented by a single data card. The format for 
these data cards is 413 ,2F1 2 .4 
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The variables read in by GENFIL are: 



NAME 


TYPE 


DESCRIPTION 


Nil 


1*4 


Mesh row nuntier of beginning of line 
or block. 


NI2 


1*4 


Mesh row number of end of line or 
block . 


NJ1 


1*4 


Mesh column number of beginning of 
line or block. 


NJ2 


1*4 


Mesh column number of end of line 
or block. 


AV 


R*4 


Normalized potential common to all 
points encompassed by Nil, NI2, 

NJ1 and NJ2. 


GV 


R*4 


EITHER: Normalized charge density 
common to all mesh points encompassed 
by Nil, NI2 , NJ1 and NJ2. 

OR: Potential description value (PDV) 
if the defined points are at a fixed 
potential . 



An arbitrary value of 99. is the PDV used to describe a point of fixed 
potential. No arithmetic operations will be carried out at this point. 
Note also that Nil < NI2, and NJ1 s NJ2, since these variables are the 
beginning and end of a DO-LOOP. In the case of a point, Nil = NI2 and 
NJ 1 = NJ2. In the case of a line of points, either Nil = NI2 or NJ1 = 
NJ2, depending on the orientation of the line. The values assigned to 
the above variables to describe the representative geometry are shown 
bel ow : 



Nil 


NI2 


NJ1 


NJ2 


AV 


GC 


Corresponding to 


1 


6 


5 


15 


1.0 


99. 


NBL0CK 1 


5 


6 


4 


4 


1.0 


99. 


NBL0CK 2 


6 


6 


3 


3 


1.0 


99. 


NBL0CK 3 


1 


4 


1 


1 


0.0 


99. 


NBL0CK 4 


1 


3 


2 


2 


0.0 


99. 


NBL0CK 5 


14 


20 


1 


1 


-0.1 


99. 


NBL0CK 6 


15 


20 


2 


2 


-0.1 


99. 


NBL0CK 7 


17 


20 


3 


6 


0.5 


0.0 


NBL0CK 8 


At the same 


time 


that the 


values 


of AV and 


GV are 


being loaded into 



appropriate arrays, the coordinates of those points with a PDV of 99. 
are being converted into plotter pen position coordinates to aid in 
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plotting the output requested by NOLINE. The pen position coordinates 
are stored in the appropriate vectors, XLIN and YLIN for later use during 
the output phase. 

Subroutine CALC1 solves the appropriate form of the finite difference 
equations, depending on the geometry variable NDIMEN, and location of the 
point in question. In addition, once the potential distribution has been 
calculated, CALC1 determines the E field, and the longitudinal and radial 
components of the E field. The value of the E field and its components 
are normalized with respect to the breakdown field strength of the 
dielectric, BRKDWN, while the values of potential are normalized with 
respect to the maximum fixed potential, GAPPOT. The derivation of the 
normalization of E is shown in Appendix A, which gives: 




The characteristic EGD geometry length ' h ' is read into the program by 
the value of the variable CONVER. Defining: 

$ALPH = r^r- 
hE b 

| E | = $ALPH | - VA | 

The value of the constant $ALPH is calculated by MAIN, and is passed to 
CALC1 as a calling argument. By normalizing E in this manner, a quick 
glance at the results will show whether or not the electric fields pro- 
duced will exceed the dielectric breakdown strength. The results of 
CALC1 are printed by subroutine WRITER. 

For the graphical output of the results, subroutine OUTPUT coordina- 
tes the various plotting subroutines required. The remainder of the 
required data cards are read in by subroutines OUTPUT under NAMELIST 
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control. The first data card is &TABL, starting in column 2. The 
variables in NAMELIST TABL are read into two vectors, TABLE and LTG, 
in the format: 

"Array name" = a set of FORTRAN constants, separated by a comma. 

The format for the elements of TABLE is TABLE( ) = ' ', where (-) 

indicates a card column available for a data item. This format prepares 
the data to be used in labeling the graphical output. The variables 



read in are: 






NAME 


TYPE 


DESCRIPTION 


TABLE (1) 


R*8 


Nozzle radius in centimeters ^ 


TABLE (2) 


R*8 


Initial charge density in C/rn 


TABLE (3) 


R*8 


Maximum fixed potential (GAPPOT) 


TABLE (4) 


R*8 


Dielectric premittivity 


TABLE (5) 


R*8 


Breakdown field strength 
of dielectric 


TABLE (6) 


R*8 


Characteristic length modeled 



The elements of the vector LTG can be used to construct a grid for the 
plotted output, label it and the plotted contours. The vector LTG is 
typed LOGICAL * 1. The variables read into the vector LTG are: 



LTG (1) 


L*1 


.TRUE. 


All exterior contour segments 
and those interior contour 








setments which represent a 
local maximum will be labeled. 






.FALSE. 


Omits labeling option. 


LTG (2) 


L*1 


.TRUE. 


Request tic marks and corres- 
ponding scale values one inch 
apart on the exterior frame 
of the contour graph. 






.FALSE. 


Omits "tic" option. 


LTG (3) 


L*1 


.TRUE. 


Request a one inch by one 
inch straight line grid to 
superimposed on the contour 
graph . 






.FALSE. 


Omits the grid option. 



Again, the last card signals the end of the NAMELIST. The format is: 

SEND 

with the symbol ampersand in column 2 of the data card. 
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APPENDIX D: MOBILITY MEASUREMENTS 



The apparatus used to investigate the mobility of the charged particles 
is shown in Figure (19). This consists of a fine Monel steel grid con- 
centric with a stainless steel cylinder. The EGD process under investiga- 
tion involves the transport of negatively charged particles. Hence if 
a voltage is applied to the grid and varies from zero to some value below 
ground, the grid will repel or pass negatively charged particles. If at 
the same time, the outer cylinder is maintained at ground, those particles 
passed by the grid will be accelerated across the gap, giving a noticeable 
variation in current flow due to the charge arriving at the outer cylinder. 
From the definition of drift velocity: 

v = k E 

and assuming that the E field is uniform: 




where 1 is the spacing between the grid and outer cylinder. But since: 

V - 1 

The mobility may be determined by: 

k . i! 

K tv 

By varying the frequency with which the grid voltage changes and ensuring 
that the waveform is a square wave, the entire range of mobility possessed 
by the charged particles can be estimated, assuming that mobility is 
proportional to the electric field between the cylinders. This is 
because at the lower frequencies, all particles will be affected, while 



44 



at the higher frequencies, or small t, inertia effects will deter all 
but the smallest or most highly charged particles from completing their 
journey. This is similar to the method of Rutherford as outlined by 
Darrow [14]. 

The measurements leading to an approximation of the mobility distri- 
bution were made at a reservoir pressure of eight psig, over a wide range 
of temperatures, 238 - 296°F. After steam flow conditions had stabilzed, 
i.e., proper pressure and constant temperature, the device was inserted 
into the EGD channel to such a position that it was axially symmetric 
with the steam jet, but not touching. This was done to keep condensation 
on the grid to a minimum. The grid was powered by a square wave genera- 
tor while the outer ring was maintained at ground. See Figure (18) for 
a schematic of the setup. Current from the ring was fed to an oscillo- 
scope for observation. It was found that the ring current varied con- 
siderably as the charged particles arrived from the Monel grid. This 
variation existed over a range of fequencies. At high temperatures, 
296°F, the range of current excitation appeared to be concentrated at 
a single point, 2.1 kHz. As the temperature of the steam was decreased, 
the range increased considerably, until at 238°F, the current excitation 
existed from an upper frequency of 5 kHz to a lower frequency of .5 kHz. 

At 296°F, the mobility corresponding to 2.1 kHz is 4.15 x 10’ 4 
m^/V-s, while at 238°F, the mobility range runs from 1 x 10~^ m^/V-s to 
1 x 10 m /V-s. The range of mobilities appears to arise from both the 
ability of large particles to move at low frequencies and the lower tem- 
perature producing a greater distribution in particle size. However, 
further work remains to be done to determine both size distribution and 
further definition of the mobility distribution. 
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Figure 1. EGD Generator Schematic. 
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Figure 2. Sample EGD Geometry and Boundary Conditions. 




Figure 3. Computer Prediction of Parallel Electrode Normalized 
Potential Distribution at Breakdown. 
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AEROSOL INJECTOR 



#316 STAINLESS STEEL ROD - 0.062" DIAMETER 




•* 0.5*- Scale - lx 



COLLECTOR PROBE 



Figure 4. Aerosol Injector and Collector Probe. 
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Figure 5. Steam Jet Schlieren Photographs 
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Figure 6. Steam Jet Profile (Axi symmetric) . 
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STEAM (FOR IONIZATION) 




o 



I 



It 



It 



52 



Figure 7. Guard Ring Schematic. 
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Figure 8. Corona Current and Convection Current Efficiency 
Ratio vs. Corona Voltage. 
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Figure 9. Potential Distribution, P = 8 psig, p = .0103 C/m 
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SLIP PARAMETER 




Figure 11. Computer Prediction of Centerline Distribution 
of Slip Parameter. 
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Figure 12. Potential Distribution, p = 0.0 C/m (Laplacian). 
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Figure 13. Electric Field Distribution, P = 14 psig, p = .0084 C/nT. 
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Figure 15. Potential Distribution, P = 14 psig. Guard Ring Voltage 
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Figure 16. Potential Distribution, P = 14 psig, Guard Ring Voltage 
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Figure 17. Potential Distribution, P = 14 psig. Guard Ring Voltage = 2 kV. 
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Figure 18. Comparison Between Laplacian |E| Field Distribution and Poisson |E| Field 
Distribution Using the Guard Ring. 
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Figure 19. EGD Mobility Measuring Schematic. 



COMPUTER PROGRAM 



CGMM0N/6L1/A(81,61 ) 

COMMON /BL2/B( 81 , 61 ) 

C0MM0N/BL3/G(81,61) 

COMMON/ BL4/NOL I NE,N1, IW, IH 
COMMON/BL5/R481) 

COMMON/ BL 6/N AT RAD, RHOZRO , PRMTVT , GAP POT, BRKDWN ,CONVER , 
1NOZPOS , NGLOSS, NROWS, NCOLS, DELU , NBLOCK , NL , NL 1 , HL 
COMMON/BLIO/NDI MEN , NDL PNT , NDLRAD 
COMMON/ 6L11/XLIN(4000) ,YLIN(4000) 

NAMEL I ST/P ARAM/ NAT RAD, RHOZRO, PRMTVT , GAPPOT , BRKDWN, 
lCONVER,NOZPOS,NGLOSS,NROWS,NCOLS,NDI ME N , D ELU , NBLOCK , 
2NL , NL1 , NOL INE , NDLPNT , N DL RAD , I W , I H 
READ ( 5 , PARAM ) 

$AL P H- GAP POT * 100. / ( BRKD WN^CONVE R ) 

$C0NST=(C0NVER/100. )**2*RH0ZR0/( PRMTVT *GAPPOT) 
WRITE(6, PARAM) 

HL=l./( (NROWS-1 )*1. ) 

WRITE(6,210) $ALPH,$CONST 
DO 10 J=1 ,NCOLS 
DO 10 1=1, NROWS 
At I ,J) =0.5 
B ( I , J ) =0. 0 
10 G( I » J ) =0.0 
19 CALL GS PRD ($CONST , 834) 

IF( NDIMEN.EQ.O) GO TO 34 
WR I T E ( 6 , 240 ) 

WRI TE ( 6 ,245 ) NROWS 

WRITE (6,230) (<J,R(J)),J=1 , NROWS) 

34 CALL G ENF I L ( £70 ) 

CALL CALC1 ($ ALPH ) 

CALL OUTPUT (870) 

GO TO 70 

»2X , 5F20 . 10 ) 

,T15,5( *R( • ,13,* )=• ,*8.5,3X) ,/) 

) 
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210 

230 

240 

245 



FORMAT ( 
FORMAT ( 
FORMAT ( 
FORMAT ( 
113,/) 

70 STOP 
END 



1 * 

, 



rT49 , * THE JET RADIUS AT STATIONS 1 THRU 
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SUBROUTINE GSPRD ( $CONS T , * ) 

COMMON /BL 3/ 0(81*61) 

C0MM0N/BL5/R (81) 

COMMON /BL6/N ATRA D, RHOZRO, PRMTVT, GAPPOT, BRK DWN , CONV ER , 
1N0ZP0S ,NGLOSS, NROWS, NCOLS,DELU,NBL OCK,NL,N LI, HL 
COMMON/ 8L10/ND I MEN , NDLPNT , NDLRAD 
Q=HL**2*$C0NST 
GLOS S = 0 . 

IF(NDIMEN) 20,10,20 

10 DO 11 I =NOZPOS , NROWS 
DO 11 J=1 , NCOL S 
IF(NGLOSS.EQ.l) GLOSS=0. 

11 G( I , Ji=( l.-GL0SS)**2*Q 
IF(NOZPOS-l) 20,19,20 

19 RETURN 1 

20 IFtNOZPOS.EQ. NDLPNT) GO TO 24 
DO 21 1=1, NDLPNT 

21 R( I )=(NDLRAD-1)*HL 
K2=NDL PNT+1 

A= ( ( NDLRAD-1 )*HL )**2 
B=( ( NATRAD-1 ) *HL )**2 
VAR A= 1 A-B ) / ( (NDLPNT-NOZPOS )*1. ) 

DO 23 I =K2 , NOZPOS 

23 R( I )={ VARA*( I-NDLPNT )+A) **.5 
IF(NDIMEN.EQ.O) RETURN 1 

24 A= ( NROWS— NOZPOS )# 1 • 

DO 35 I=NOZPOS, NROWS 
B = ( I — NOZPOS) *1 . 

X=B/A 

35 R ( I ) = . 0445+. 0043*X+. C533*X**2 
B=R( NOZPOS) 

DO 50 I=NOZPOS, NROWS 
DO 40 J= 1 , NCOL S 
K1=J-1 

I F ( R ( 1 )-Kl*HL) 43,42,42 

42 G ( I , J>=(B/R(I)-GL0SS)**2*Q 
GO TO 40 

43 IF(R( I )-HL*(Kl-.5) ) 50,44,44 

44 G ( 1 , J)=CB/R(I)-GLOSS )**2*Q 
GO TO 50 

40 CONTINUE 

50 CONTINUE 

IF ( NOZPOS. EQ. NDLPNT) RETURN 

K1=NATRAD-1 

K3=2*N ATRAD 

DO 90 I=K2, NOZPOS 

DO 80 J=NDLRAD , K3 

IF ( NGLOSS. EQ. 1 ) GL0SS=0. 

NRAD=J-NDLRAD 
IF(Kl-J) 60,51,51 

51 IF( R( I )— NRAD*HL ) 53,52,52 

52 G( I, J)=(B/R( I)-GLOSS )**2*Q 
GO TO 60 

53 I F ( R ( I )— HL* ( NRAD—. 5 ) ) 80,54,54 

54 G(I,J)=(B/R(I)-GL0SS)**2*Q 
GO TO 60 

60 NRAD=2*NDLRAD-J 
IF ( NRAD ) 80,80,65 
65 G( I ,NRAD)=(B/R(I )-GL0SS)**2*Q 
80 CONTINUE 
90 CONTINUE 
RETURN 
END 
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SUBROUTINE GENFILt*) 

COMMON/BL 1/A(81,61) 

COMMON/ BL2/B (81 ,61 ) 

C0MM0N/8L3/G ( 81,61) 

COMMON/ B L4/N0L INE,N1,IW,IH 

COMMON /B L6 /N AT RAD, RHOZRO, PRMTVT, GAPPOT, BRK OWN, CON VER, 
1NOZPOS ,NGLOSS,NROWS,NCOLS,DELU,NBLOCK,NL,NLl ,HL 
COMMON/BL ll/XL IN (4000) ,YLIN(4000 ) 

XVAL=( IW*1.)/(NC0LS*1.-1.) 

YV AL= { IH*1.)/(NR0WS*1.-1.) 

N1 =0 

DO 20 1=1 , NBLOCK 

READ (5, 200) N 1 1 , N 1 2 , N J 1 ,N J2 , P V , GV 
DO 10 I1=NI1,NI2 
DO 10 J1=NJ1,NJ2 
A ( I 1, J1)=PV 
B ( 11 , J1 )=PV 
G{ 11 ,Ji)=GV 
I F ( GV— 98 .98 ) 10,5,5 
5 N1 =N1 + 1 

XLIN(N1)=XVAL 5|; (J1-1) 

YLIN(N1)=YVAL*{ I 1- 1 ) 

10 CONTINUE 

20 CONTINUE 
IF(NOLINE) 22,21,22 

21 CALL OUTPUT ( 870 ) 

70 RETURN 1 

22 RETURN 

200 FORMAT (41 3,2F12.3) 

END 
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SUBROUTINE CALCltSALPH) 

COMMON /BL1/A181, 61) 

C0MM0N/BL2/B(81,61) 

C0MM0N/BL3/G( 81,61) 

COMMON/BL5/RI81) 

COMMON /BL6/ NAT RAD, RHQZRO, PRMTVT , GAPPOT , BRK DWN , CONV ER , 
1N0ZP0S,NGL0SS,NR0WS,NC0LS,DELU,NBL0CK,NL,NL1 ,HL 
COMMON/ BL 10/ND I MEN, NDLPNT,NDLR AD 
I TCNT= 1 
ITERS= 1 

5 DO 20 J = 1 , NCOL S 
J1=J+1 
J2=J-1 
R1=J2*HL 
IF(R1) 7,7,8 

8 ARG 1= . 5*HL /R 1 

7 DO 15 1=1, NROWS 
11=1+1 
12=1-1 

IF(G(I,J).GT.98..AND.G(I,J).LT.100.) GO TO 15 
IFt J.EQ.l.OR.J.EQ.NCCLS) GO TO 11 
IFU.EQ.l.OR.I.EQ. NROWS) GO TO 16 
IF ( NDIMEN) 10,10,9 

9 At I ,J)=B( I ,J)+1»5*(.25*(G( I jJJ+ARGl^tAt I , J 1 ) -A ( I,J2) ) 
1+A( I ,J1)+A(I,J2)+A( I1,J)+A(I2, J) )-B( I,J) ) 

GO TO 15 

10 At I , J) =B{ I , J ) + l. 5*t . 25* tGt I , J) + A{ I , Jl) +A ( I , J2)+At 1 1 , J) 
1 +A ( 12, J) ) -81 I, J) ) 

GO TO 15 

11 IF{ J .EQ.NCOLS) GO TO 12 

A{I,J)=B{ I,1)+1.5*(1.333*A( I,2)-.333*A{ I,3)-B{ 1,1) ) 

GO TO 15 

12 J3= j-2 

At I, J) = B( I ,J )+1.5*( 1.333*A( I , J2 )-. 333* At I , J 3 ) -B ( I , J ) ) 
GO TO 15 

16 IF ( I .EQ. NROWS) GO TO 18 

Ail, J)=B(l,J)+1.5*ll.333*At2,J )-.333*A ( 3. J )-Bt 1, J ) ) 

GO TO 15 
18 13=1-2 

At I , J) = Bt I , J)+1.5*( 1 .333*A { 12, J )-. 333*A( I3,J)-B( I,J) ) 
15 CONTINUE 
20 CONTINUE 
DELMAX=0. 

DO 24 J = 1 ,NCOL S 
DO 24 1=1, NROWS 
DEL=A( I , J ) — B t I , J ) 

ADEL=ABS { DEL ) 

IF ( ADEL.GT.DELMAX) DFLMAX=ADEL 
24 CONTINUE 

IF ( DELU-DELMAX) 26,40,40 
26 DO 34 J = 1 , NCOL S 
DO 34 1=1, NROWS 

34 Bt I, J)=A{ I ,J) 

IF C ITCNT-100) 36,35,35 

35 WRITE(6,210) ITERS, DELMAX 
ITCNT= 1 
ITERS=ITERS+1 

GO TO 5 

36 ITERS=ITERS+1 
I T CNT = I TCNT + 1 
GO TO 5 

40 WRITE{6, 320 ) ITERS 
WR I T E 1 6 , 230 ) 

WR I TE t 6 , 232) 

CALL WRITER! A, NROWS , NCOL S) 

ARG2 = . 5/HL*$ ALPH 
DO 45 J= 1 , NCOL S 
J1=J+1 
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J2= J— 1 
J3=J + 2 
J4= J-2 

DO 45 1=1, NROWS 
11=1+1 
12=1-1 
13=1+2 
14=1-2 

IF! I .EQ.l. OR. I .EQ.NROWS) GO TO 41 
AZ=A < I 2 , J ) —A ( II, J) 

GO TO 42 

41 IFU.EQ.l) AZ=A( 3,J ) +3 . *A { 1 , J )-4 . *A ( 2,J) 

I F { I .EQ.NROWS) AZ=4.*A( 12, J)-A( I 4 , J ) -3 . *A ( I , J ) 

42 IFU.EQ.l. OR. J.EQ.NCOLS) GO TO 43 
AR=A { I , J 2 ) —A { I , J 1 ) 

GO TO 44 

43 AR= 0.0 

44 B( I , J ) =ARG2*AZ 
G ( I , J ) =ARG2* AR 

45 CONTINUE 

WR ITE(6,230) 

WRI TE ( 6 , 237 ) 

CALL WRITERtB, NROWS, NCOLS) 

WRITE (6,230) 

WRITE! 6,239) 

CALL WRITER!G, NROWS, NCOLS) 

EMAX = 0 . 0 

DO 55 J=1,NC0LS 

DO 55 1=1, NROWS 

SUM= B { I , J ) **2 + G ( I, J)**2 

B { I , J ) =SQRT( SUM) 

IF! B( I , J ) .GT.EMAX) EMAX=B(I,J) 

55 CONTINUE 

WRITE { 6,230) 

WRITE! 6,234) 

WRITE!6,235) EMAX 

CALL WRITER!B, NROWS, NCOLS) 

IP(NDIMEN) 58,58,109 

109 DO 49 1=1, NROWS 

I F ( R { I ) ) 110,110,120 

110 DO 113 J=l, NCOLS 
113 G! I , J) =1.+1.*{ J-l) 

GO TO 49 
120 XS L P = 1 ./R! I ) 

DO 49 J=l, NCOLS 
G ( I ,J )=XSLP*HL*!J-1) 

49 CONTINUE 

210 FORMAT!' • , T5 , • A FT ER • , I 5 , 1 X , ' I T ERAT I ONS — DELMAX= ' , 

1 FI 1.5,/) 

230 FORMAT Cl') 

232 FORMAT!' »,T49,'THE MATRIX OF POTENTIALS',/) 

234 FORMAT!’ ',T51,'THE |E| FIELD MATRIX',/) 

235 FORMAT!' • ,T42 ,' MAXI MUM FIELD STRENGTH IS',F12.4,/) 
237 FORMAT!' ',T46, 'AXIAL COMPONENTS OF |E| FIELD',/) 
239 FORMAT!' • ,T47 , • RADI AL COMPONENTS OF |E| FIELD',/) 
320 FORMAT! • ' , ' THE POTENTIAL ITERATION SCHEME HAS', 

1T37, 'CONVERGED IN', 15,' ITERATIONS',/) 

58 RETURN 
END 
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SUBROUTINE WR I TER ( Z , NROWS , NCOL S ) 
DIMENSION Z(81,61) 

J1 = I 
J 2= 1 0 

10 DO 15 1=1, NROWS 

IF (J2.GE. NCOLS) J2=NC0LS 
15 WR I TE I 6 , 220) I ,(Z( I , J ) ,J=J1,J2) 
IF (J2 . EQ. NCOLS ) RETURN 
J1=J1+10 
J2=J 2+10 
WRITE (6,230) 

GO TO 10 

220 FORMAT ( ' ',13,10^12.4) 

230 FORMAT ( ' 1' ) 

END 



SUBROUTINE CLVL ( CM , CLM , NUML ) 

C0MM0N/BL6/NATRAD, RHCZ RO , P RMT VT , GAP POT , BRKDWN , CONV ER , 
1N0ZP0S , NGLOSS, NROWS, NCOL S, DE LU, NBLOCK, NL, NL1 ,HL 
DIMENSION CM (81, 61 ),CLM(NUML) 

CMIN=CM( 1,1) 

CMAX=CMIN 
DO 5 J = 1 , NCOLS 
DO 5 1=1, NROWS 

IF (CM( I , J ) .LT.CMIN ) CMIN=CM(I , J) 

IF (CM ( I ,J) .GT.CMAX) CMAX=CM(I,J) 

5 CONTINUE 

C NOW DETERMINE CONTOUP LEVELS TO BE PLOTTED. 

J=NUML 
I=NUML-1 
ANL = 1*1. 

PL I NT= ( CMAX-CM I N ) / ANL 
C NOW FILL THE CONTOUR LEVEL VECTOR 

DC 6 1=1, J 

6 CLM(I ) =CMIN+(I-1)*PLINT 
RETURN 

END 



SUBROUTINE F LOP ( Z, NROWS , NCOLS ) 
DIMENSION Z(81,61) 

1 1 VRT = NROW S/ 2 
DO 3 1=1,1 IVRT 
M=NROW S — < I — 1 ) 

DO 3 J=l, NCOLS 
SAVE=Z ( I , J) 

Z( I, J)=Z(M,J) 

3 Z(M,J)=SAVE 
RETURN 
END 
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SUBROUTINE OUTPUT! 
REAL*8 TI TIE ( 30) / ' 



1 ' 


. , . , 


, 


* , ' NOZZLE 


• 


2 ' 




, • PRMTVT 


, t 


1 


3 • 


' , ' RHOZRO • 


• 


• t 


1 


4* GAPPOT 




i 


• , • BRKDWN 


1 


5' 


' ,' ' 


, * LENGTH 


i t 


• 


6 ' 


* , ' C . M . BOH' 


, 'LEY 


*,' POTENTI 


« 


7 * AL MAP 




, ' BOX 


89'/ 




REAL*8 


T I TLE A ( 2 ) / ' IE | 


FI EL * , ' 


D MAP '/ 




REAL *8 


TITLEB ( 2 ) / * 


GEOMET* , 


•RY PLOT ' / 




REAL*8 


TABLE (6 ) 









LOG I C A L* 1 LTG ( 3 ) 

CO MM ON/ B LI/ A ( 81, 61 ) 

C0MM0N/BL2/B(81 ,61 ) 

COMMON /BL3/G ( 81,61) 

C0MM0N/BL4/N0LINE,N1, IW, IH 
COMMON /BL5/R(81 ) 

COMMON /BL6/N ATRAD, RHCZRO , PRMTVT ,GAPPOT , BRKDWN ,CONVER , 
1N0ZP0S , NGLOSS, NROWS , NCOLS , DELU,NBL0CK,NL,NL1,HL 
COMMON /BL7/CL(30) 

COMMON /BL8/C LI ( 30) 

COMMON /BL9 /CL 2 ( 1 ) 

NAMELI ST/TABL/ TABLE, LTG 
RE AD ( 5 , T ABL ) 

FIRST FIGURE CONTOUR LEVELS TO BE PLOTTED FOR MATRIX A 
CALL CLVL ( A, CL ,NL) 

NOW FIGURE CONTOUR LEVELS FOR |E| FIELD. 

CALL C LVL ( B , C L 1 , ML 1 ) 

THE CONTOUR LEVEL VALUE FOR THE CLOUD BOUNDARY IS 1.0 
CL 2 ( 11-1.0 
WRI TE ( 6 ,210) 

DO 3 1=1, NL 

l WRITE(6,215) I, CL (I), I, CL 1(1) 

CALL FLOP (A, NROWS, NCOLS) 

CALL FLOP (B, NROWS , NCOLS) 

CALL FLOP (G, NROWS, NCOLS ) 

TI TLE ( 8 ) =TABLE ( 1 ) 

T I TL E ( 11)=TABLE( 4) 

TITLE ( 14 ) =TABL E ( 2 ) 

TITLE( 17) =TABLE(3) 

T I TL E ( 20 ) =TABL E ( 5) 

TITLE(23 )=TABLE(6) 

IF(NOLINE) 7,6,7 
» TITLE! 27)=TITLEB( 1 ) 

TITLE(28)=TITLEB(2) 

’ CALL CONTURtA, NROWS, NCOLS, 81, CL, NL , TITLE, LTG, G,CL2, 

11 , £ 20 ) 



NOW PLOT |E 
TITLE! 27)=T 
T IT L E ( 28 ) =T 



FIELD CONTOURS. 

TLE A! 1 ) 

[TLE A! 2 ) 

CALL C ONTUR(B, NROWS, NCOLS, 81 ,CL1,NL1,T ITLE , LTG, G, 
1CL 2, 1, 820) 

210 FORMAT {• *,T26,'THE POTENTIAL LEVELS PLOTTED ', T80 . 

1 ' THE | E | " 

215 FORMAT (• 

1F8.4,/) 

20 IF(NOLINE) 

10 RETURN 1 
12 RETURN 
END 



FIELD 

' , T33 , 



LEVELS PLOTTED* ,//) 

•CL! • ,1 2, * ) =' ,F8.4,T87,'CL1( • , 12, ' )=' , 



12 , 10,12 
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SUBROUTINE CONTUR ( AM , M ,N , MX , CL , NL, T I TL E, LTG , G , CL 2 , NL2, 
1 * ) 

REAL *8 TITLE* 1 ) 

REALMS WIDTH/' WIDTH' /, HEIGHT/' HE IG HT ' / , WH I CH 
DIMENSION AM ( MX » 1 ) ,CL( 1) 

DIMENSION G( MX , 1 ) , CL2( 1 ) 

DIMENSION RECI900), X(1800), Y(1800) 

DIMENSION I P T ( 3 , 3 ) , I NX( 8 ) , I NY ( 8 ) 

COMMON /DAYHOF/ MT , NT , N I , I X , I Y , I DX , I DY , I SS , I T , I V , NP , 
1NQ, JT,PY,REC,CV, IPT, I NX, IN Y , DL , R A , THE 
COMMON /INTFAC/ X,Y 
DIMENSION DITSX(5),DITSY(5 ) 

LOG I CAL* 1 LTG( 1 ) , MINUS ,LABL 
COMMON /TAB L/ T A BC ( 20 , 6 ) , JC 

COMMON/D I T S/XM I N,Y MIN, SLOP EX, S LOPE Y, DI TSDX , D I T SDY , 
IIDIR,LABL , MINUS 
COMMON /BL4/N0L I NE,NI, IW, IH 
COMMON/BL10/NDIMEN,NDLPNT, NDLRAD 
C0MM0N/BL13/NSCN 
JC = 0 

LABL=LTG ( 1 ) 

C CHECK IW PARAMETER 

WH ICH=W I DTH 
IF(IW) 1,1,2 

1 WR I T E ( 6 » 60 ) WHICH 

60 FORMAT ( 'O' ,T7, A8, 'OF CONTOUR GRAPH ILLEGAL.') 

71 WR I TE ( 6 , 64 ) 

64 FORMAT* 'O' ,T7, 'NO GRAPH WILL BE PRODUCED.') 

RETURN 

C CHECK IF IW IS TOO WIDE 

2 IF ( IW-9) 3,3,40 

40 WR ITE(6»61) 

61 FORMAT ( ' O' ,T7, • I W PARAMETER GREATER THAN 9. CONTUR 
1 WILL SET I W = 9 . ' ) 

I W=9 

C NOW CHECK IH PARAMETER 

3 IF(IH) 4,4,5 

4 WhICH=HEIGHT 
GO TO 1 

5 DITSDX= (N-1.0 ) / I W 
DITSDY=(— 1 .0+M ) / I H 
XM I N = 1 . 0 

YM IN=-M 

SL0PEX=1.0/DITSDX 
SL0PEY=1.0/DITSDY 
DI TS X ( 1 ) = 1 .0 
DI TSX( 4) =1.0 
DITSXI 5 ) = 1 • 0 
DITSX(2)=N 
DI TS X( 3 ) =N 
DITSYt 1 ) =- 1 . 0 
DI TS Y ( 2 ) =—1.0 
DITSYl 5) =-1.0 
DI TS Y ( 3 ) =— M 
D IT SY ( 4 ) =— M 
DO 2011 1=1,5 

DITSX( I )=SLOPEX*(DITSX( I )-XMIN) 

2011 DITSYl I )=SLOPEY*(DITSY (I J-YMIN) 

STARTP=(9.0-IW)/2.0 
CALL PLOTS 

CALL PL0T(STARTP,0.0,-3) 

CALL LINE(DITSX,DITSY»5»1,1) 

DITSX( 1)=DITSX( 1 )— .5 
DITSX(5)=DITSX(1) 

DITSX(4)=DITSX(4)-.5 
DITSX(2)=DITSX(2)+.5 
DITSX(3)=DITSX(3)+.5 
DI TS Y ( 1 ) =D I TSY ( 1 )+.5 



72 



DITSY (5 )=DITSY { 1 ) 
DITSY(2)=DITSY(2)+.5 
DITSY( 3)=DITSY(3)-.5 
DITSY(4)=DITSY (4)-. 5 
CALL LINE(DITSX,DITSY,5,1,1) 
SL0PEX=1.0/DITSDX 
SL0PEY=1.0/0ITSDY 
IENDX=SL0PEX*N+1 
IENDY=SL0PEY*M+1 
I F ( . NOT. LTG ( 2 ) ) GO TO 34 
C DRAW TIC MARKS QN OUTER FRAME 

C START ON LEFT EDGE GOING DOWNWARD 

I F LAG = 0 
Z I NGX=-. 1 
ZI NGY = 0 . 0 
ZX=0.0 
ZY=- 1.0 
CX = D IT SX ( 1 ) 

CY=D ITSY(l) — .5 
I END=I ENDY 
2222 IFLAG= IFLAG+1 

DO 2022 1 = 1, I END 
CALL PLOT (CX,CY,3) 

COORDX=CX+ZINGX 
COORDY=CY+ZINGY 
CALL PLOT (COORDX,COORDY,2) 

CX=CX+ZX 
2022 CY=CY+ZY 

GO TO ( 21,22,23,24) , IFLAG 
C NOW DO THE RIGHT EDGE GOING DOWNWARD 

21 Z I NG X= . 1 

CX = DIT SX ( 2 ) 

CY=DITSY(2)-.5 
GO TO 2222 

C NOW DO TOP EDGE 

22 ZI NGX=0. 0 
Z I NGY= • 1 
ZX = 1.0 
ZY=0. 0 

CX=D I T SX { 1 ) + . 5 
CY = D I TS Y ( 1 ) 

I END=I ENDX 
GO TO 2222 

C NOW DO THE BOTTOM EDGE 

23 Z I NGY=— . 1 
CX=D IT SX ( 4 ) + . 5 
CY=D I T SY ( 4 ) 

Z I NGY=-. 1 
GO TO 2222 

C NOW LABEL TIC MARKS 

C DO X-DIRECTION FIRST, TOP EDGE 

C POSITION PEN 

24 DE LT AX =DI T SDX 
I FLAG=0 

ZX = 1 .0 
ZY=0. 0 

CX=DITSX( l)+.35 
CY=DITSY(1)+.12 
3033 I F LAG= I FLAG+ 1 
XZERO= 1.0 
DO 3333 1 = 1, I END 

CALL NUMBER (CX,CY, . 14 , XZERO , 0 . 0 , 1 ) 

CX=CX+ZX 

CY=CY+ZY 

3333 XZ ERO = XZERO+DE LTAX 

GO TO (31,32,33,34) , IFLAG 

C LABEL BOTTOM EDGE TIC MARKS 

31 CX=DITSX{4)+.35 

CY=DITSY(4)-.19 
GO TO 3033 

C LABEL LEFT EDGE OF TIC MARKS 

32 CX = DITSX(4)-.4 
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CY=D ITSY{4 ) + .46 
DELTA X=D I TSDY 
I END= I ENDY 
ZX= 0.0 
ZY= 1.0 
GO TO 3033 

C NOW LABEL RIGHT EDGE TIC MARKS 

33 CX=D ITSXI 3)+. 12 
CY=DITSY(3)+.46 
GO TO 3033 

C CHECK IF GRID DESIRED 

34 CALL RESTOF(LTC, IENDX, I END Y , NL , AM , M , N , MX , C L , NL2 
1STARTP,TITLE,DITSX,DITSY) 

IF(NOLINE) 2502 t 2505, 2502 
2502 RETURN 
2505 RETURN 1 
END 



r* 

t & 



SUBROUTINE RESTOF(LTG» I END X , I END Y, NL » A M , M , N , MX , CL 
2,TITLE ,DITSX,DITSY) 

REAL*8 TITLE(l) 

DIMENSION AM ( MX » 1 ) ,CL( 1) 

DIMENSION G ( MX , 1 ) » CL2 ( 1 ) 

DIMENSION DITSXI 5) tDITSYI 5) 

L0GICAL*1 LTG(1),MINUS»LABL 
COMMON /TAB L/T ABC ( 2 0 » 6 ) *JC 

COMMON /D I T S/XM I N, YMI N, SLOP EX, SLOPE Y t DI TSDX , DI T SDY 
1 IDIR ,LABL, MINUS 
COMMON /BL4/N0LI NE, N1 , I W, IH 
COMMON/ BL 10/ND I M EN , NDL PNT , NDLR AD 
C0MM0N/BL1 1/XL IN (4000) ,YLIN(4000) 

COMMON /8L 13/NSC N 
I F ( .NOT .LTG(3) ) GO TO 35 
C DRAW INCH BY INCH GRID 

I END=I ENDX-2 
C POSITION PEN 

I F LAG = 0 

CX=D IT SX( 1 ) + • 5 
CY=DITSY(l)-.5 
COORDX = 0 . 0 
COORD Y=- I H 
DX = 1 .0 
DY=0. 0 

4044 DO 4444 I=1,IEND 
CX=CX+DX 
CY=C Y+DY 

CALL PLOT ( CX , CY , 3 ) 

ZX=CX+COORDX 
ZY=CY+ COORD Y 
4444 CALL PLOT ( ZX , ZY , 2 ) 

I F ( I FLAG ) 35,42,35 



CL2 , 



NL2 , 
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42 IFLAG= 1 

I END= I ENDY-2 
CY=D ITSV(4}+.5 
CX=DITSX(4)+.5 
COORDX = I to 
COORDY =0.0 
CX=0.0 
DY = 1 .0 
GO TO 4044 
35 CONTINUE 

IF(NOLINE) 11,2500,2500 
2500 CALL LINE(XLIN,YLIN,N1,1 , — 6 ) 

IF(NOLINE) 11,779,11 
11 DO 20 1=1, NL 
20 CALL SCAN! AM , M , N , MX , CL ( I ) ) 

NSCN= 1 

IF ( NDI MEN ) 2002,2003,2002 

2002 DO 2004 1=1, NL2 

2004 CALL SCAN(G,M,N,MX,CL2(I)) 

2003 I F ( .NOT.LABL) GO TO 778 
IF(JC.EQ.O) GO TO 778 
DO 777 1=1, JC 
CC0RDX=TA3C( I ,4) 

COORDY=TABC( I ,5) 

CL EV = T ABC ( 1,6) 

CALL NUMBER! COORDX, CCORDY , . 07 , CL EV , 0 .0 , 3 ) 

777 CONTINUE 

778 CALL SYMeOL(-STARTP, IH+1.0, .21, TITLF(25) ,0.0,48) 
CALL SYMBOL(-STARTP, IH+1 .5,.21,TITLE!19),0.0,48) 
CALL SYMBOL!- START P, IH+2.0 , .21 .TITLE ( 13 ) ,0.0,48) 
CALL SYMBOL! -ST ART P, IH+2.5, .21,TITLE(7),0.0,48) 
CALL SYMBOL(-STARTP, IH+3.0 , .21 .TITLE! 1) ,0.0,48) 
CALL PLOT (-STARTP, IH+6.5,-3) 

CALL PLOTE 
RETURN 

779 CALL SYMBOLt-STARTP, IH+1.0,. 21, TITLE!25) ,0.0,48) 
CALL P LOT ( —START P , IH+4 .5,-3) 

CALL PLOTE 

RETURN 

END 
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SUBROUTINE SC AN ( AM, M, N ,MX. CL ) 

DIMENSION AM ( MX , 1 ),REC<900), X(1800), Y(1800) 

DIMENSION IPT(3,3),INX(8),INY(8) 

COMMON /DAYHOF/ MT , NT, N I , I X , I Y , I DX, IDY , I SS , I T, I V , NP , NQ 

I PY,REC,CV, IPT,INX,INY,DL,RA,THE 
COMMON / INTFAC/ X,Y 

LOGICAL *1 LABL, MINUS 

COMMON /D I T S/XM I N,YMIN,S LOPE X,S LOPE Y , DI TS DX , D ITSDY , 

I I DIR , LABL , MINUS 
COMMON /BL 13/ NS CN 
D= 0. 

R=l. 

TH = 1 .570796 

NP = 0 

DL = D 

RA = R 

THE=TH 

MT =N 

NT = M 

CV = CL 

IF (NSCN.EQ.l) IZW=0 
I F ( I ZW-1 20631 ) 1,3,1 

1 I PT ( 1 , 1 ) =8 
IPT(1,2)=1 
I PT ( 1 , 3 ) =2 
I PT ( 2 , 1 ) =7 
I PT ( 2 , 3 ) =3 
I P T ( 3 , 1 }=6 
I PT ( 3 , 2 ) =5 
IPT( 3, 3 ) =4 
I NX ( 1 ) =— 1 

I NX ( 2 ) =— 1 
INX( 3 ) =0 
I MX ( 4 ) = 1 
I NX ( 5 ) =1 
I N X ( 6 ) = 1 
I NX ( 7 ) =0 
I NX ( 8) =-l 
IN Y ( 1 ) =0 
I MY ( 2 ) = 1 
INYO) =+l 
I NY ( 4 ) =+ 1 
I NY (5 )=0 
INY(6)=-1 
IN Y ( 7 ) =— 1 
I NY ( 8 ) =-l 
I Z W= 12 0631 
3 XT = MT 

DO 58 J = 1 ,900 
58 REC ( J ) =0 
ISS = 0 

2 MT1 =MT-1 
I D I R = 1 

DO 110 1=1, MT1 
IF(AMt 1,1)— CV) 55,110,110 
55 IF(AM( 1,1 + 1)— CV) 110,57,57 
57 IX=I+1 
I Y =1 
I D X=— 1 
I DY = 0 

CALL TRACE (AM, MX) 

110 CONTINUE 
NT1=NT— 1 
I D I R =2 

DO 2 0 I = 1 , NT 1 
I F ( AM ( I , MT )— CV ) 15,20,20 

15 IF C AM( I + 1,MT)-CV) 20,17,17 
17 I X=MT 
I Y = I + 1 
I DX = 0 
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IDY=- 1 

CALL TRACE (AM, MX) 

20 CONTINUE 
I D I R = 3 

22 DO 30 1=1, MT1 
MT 2 = MT + 1—1 

IF (AMI NT ,MT2 ) — CV ) 25,30,30 
25 IF(AM(NT,MT2-1 )-CV) 30,27,27 
27 IX=MT2— 1 
IY=NT 
IDX = 1 
IDY=0 

CALL TRACE (AM, MX) 

30 CONTINUE 
ID IR=4 

DO 40 I = 1 , NT 1 
NT2=NT+1— I 

IF ( AM( NT2, 1)-CV ) 35,40,40 
35 IF(AM(NT2-L,1)-CV) 40,37,37 
37 IX=1 

I Y = NT 2—1 
I DX=0 
I DY= 1 

CALL TRACE (AM, MX) 

40 CONTINUE 
I D I R = 5 
I S S= 1 
NT1 =NT— 1 
MT 1=MT— 1 
DO 1C J = 2 , NT 1 
DO 10 1=1, MT1 
I F ( AM ( J, I )-CV) 5,10,10 
5 IF C AM f J, I + 1)-CV) 10,7,7 
7 C0M=1C0*( 1+1 )+J 
IF (NP) 12,11,12 
12 DO 9 I D= 1 , NP 

IF (REC( ID)-COM) 9,10,9 
9 CONTINUE 
11 IX= 1+1 
I Y = J 
IDX=- 1 
I DY = 0 

CALL TRACE (AM, MX) 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE TRACE (AM, MY) 

DIMENSION AM (MY , 1 ),PEC(900), X(1800), Y(1800) 

DIMENSION I PT ( 3 , 3 ) , INX( 8 ) , INY ( 8 ) 

COMMON /DAYHOF/ MT , NT , N I , I X , I Y , I DX , I DY , I SS , IT , I\/ , NP , 
IN »JT,PY»REC,CV,IPT,INX»INY,DL»RA»THE 
COMMON / INTFAC/ X,Y 
PY= 0. 0 

RC= COS ( THE ) *RA 
RS= SIN ( THE ) *R A 
501 J T = 0 
N=0 
IXO=IX 
I Y0= I Y 
IS X= I DX+2 
ISY=IDY+2 



I S= I PT ( I SX , I SY ) 

JT B= 0 
ISO=IS 

I F ( ISO-8) 18,18,17 

17 I S 0= I S 0- 8 

18 IT=0 

5 CONTINUE 

CALL CALC (AM, MY) 

NZ = N 
N=NZ 

IF (IT+JT-1) 49 ,49,47 
47 XS=X ( N- 1 ) 

YS=Y ( N-l ) 

X(N-1)=X(N) 

Y ( N- 1 ) =Y ( N ) 

X(N) =X S 

Y ( N ) =Y S 
49 IS=IS+1 

JT = I T 

9 IF ( I S — 9 ) 8,7,7 

7 IS=I S-8 

8 IDX=INX(IS) 

IDY=INY( IS) 

IX2=IX+IDX 
IY2=I Y+IDY 
JTB= JTB+1 

IF (JTB-1799) 51,51,308 
308 PRINT 1 03 , CV , X ( N ) ,Y(N) 

103 FORMAT ( 1H0,23HA CONTOUR LINE AT LEVEL,E12.5, 

121 H WAS TERMINATED AT X=,E12.5,3H Y=,E12.5/ 

2 48H BECAUSE IT CONTAINED MORE THAN 1799 PLOT POINTS 
RETURN 
51 CONTINUE 

IF (ISS) 10,10,20 

20 IF(IX-IXO) 12,21,12 

21 IF(IY-IYO) 12,22,12 

22 IF(IS-ISO) 12,23,1 2 

23 CONTINUE 

CALL CALC (AM, MY) 

GO TO 73 

10 I F ( I X2 ) 13,50,13 

13 IF (IX2-MT) 19,19,50 

19 IF ( IY2) 11, 50,11 

11 IF (IY2-NT) 12,12,50 

12 I F ( C V- AM ( I Y2 , I X2 ) J 206,206,5 

206 IF { IDX**2+IDY**2-1) 213,6,213 

213 DCP=(AM(IY,IX)+AM( IY, I X2 ) +AM ( I Y2 , 1 X ) +AM( I Y2 , 1 X2) ) /4. 0 
IF (DCP-CV) 5,217,217 

217 IF (INX(IS-l)J 214,215,214 

214 IX= I X+ I DX 
I DX=- I DX 
PY=2 . 0 

CALL CALC (AM, MY) 
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I X= I X+IDX 
GO TO 6 
215 I Y= I Y+ I DY 
I D Y=— I DY 
PY = 2 . 0 

CALL CALC (AM,m 
IY=I Y+IDY 

6 IF( AM( IY, IX-1 )-CV) 306,16,16 

306 NP=NP+1 

REC( MP )=100*IX+I Y 
16 I$= I S+5 
IX=I X2 
I Y= I Y2 
GO TO 9 
50 XT =MT 

I F ( AM ( IY, IX-D-CV) 307,73,73 

307 NP=NP+1 

REC ( NP ) = 1 00* I X+ I Y 

73 DO 74 I = 1 , N 

X ( I ) =X ( I ) + RC-Y ( I ) 

74 Y { I ) =RS*Y( I ) 

CALL PLOTT l N , C V ) 

RETURN 

END 
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SUBROUTINE CALCIAM,MY) 

DIMENSION AM(MY,1 ),REC1900), XI1800), YI1800) 

DIMENSION IPT (3*3) , I NX I 8 ) , INY( 8) 

COMMON /DAYHOF/ MT , NT , N I , I X , I Y , I DX , I DY , I SS , IT , W , NP , N, 
1JT,PY,REC,CV,IPT,INX,INY,DL,RA,THE 
COMMON /INTFAC/ X,Y 
I T = 0 
N=N+ 1 

IF <IDX**2 + IDY**2 -1) 20,1,20 

1 IF (ICX) 10,2,10 

2 X { N ) = I X 
Z= I Y 

I Y2= I Y+ I DY 
DY= I DY 

41 Y ( N ) = ( { AM ( IY, IX ) — C V ) / ( AMI I Y, IX ) — AM { IY2, IX) ) )*DY+Z 
RETURN 
10 YI N ) = I Y 
W=IX 
DX=I DX 
I X2= I X + 1 DX 

44 X(N) = ( (AMI IY, I X ) — C V ) / { AMI IY, IX) -AM I IY, 1X2) ) )*DX+W 
RETURN 

20 I X2= I X + I DX 
I Y2 = I Y + 1 DY 
H=IX 

Z= I Y 
DX= I DX 
DY= I DY 

DCP= I AMI IY, IX)+AM( IY, IX2J + AMI IY2,I XJ+AMI IY2, 1X2) ) /4. 0 
IF IPY-2.0) 24,21,24 

24 IF (DCP-CV) 21,21,25 

21 AL = AM I IY, IX)-DCP 

23 V=. 5*1 AL+DCP -CVJ/AL 

27 XIN)=V*DX+W 
YIN) =V*DY+Z 
PY=0. 0 
RETURN 

25 IT=1 

AL = AM I I Y2 ,1X2 )-DCP 
33 V= . 5* I AL+DCP— CV ) /AL 

28 X ( N ) =- V*DX+W + DX 
Y(N)=-V*DY+Z + DY 
YIN) =- V*DY +Z + DY 
RETURN 

END 
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SUBROUTINE PLOTT(NP,CV) 

COMMON /INTFAC/X( 1800) ,Y( 1800) 

LOG I CAL* 1 MINUS, LABL 
COMMON /TAB L/ TABC ( 20 , 6 ) , JC 

COMMON /D I T S/XM I N,YMIN, SLOP EX, SLOPE Y,DI TSDX , D I T SD Y , 

1 IDIR,LABL, MINUS 

C SCALE POINTS FOR PLOT ROUTINE 

DO 100 1=1, NP 
X< I )=SLOPEX*(X ( I J-XMIN) 

100 Y( I )=SLOPEY*(-Y( I I-YMI N) 

CALL LINE(X,Y,MP,1,1) 

IF ( . NOT. LABL) RETURN 
D I R = 0 . 0 

GO TO (1,2,3,4,61, ICIR 

1 D I R = 9 0 . 

2 COCRDX=X ( 1 ) 

CGC'RDY =Y ( 1 ) 

5 CALL NUMBER ( COORDX , CCORDY , . 07 , CV , 0 I R , 3 ) 

RETURN 

C MOVE PEN DOWN ONE HALF INCH 

3 DIR=90. 

COORDX=X ( 1 ) 

COORDY =Y ( 1 ) - . 3 
GO TO 5 

C MOVE PEN TO THE LEFT 

4 COORDX = X ( 1 ) - . 3 
COORDY =Y ( 1 } 

GO TO 5 

C SEARCH FOR XM AX , XM IN , YMAX , YM I N , AND SAVE YM I NX 

6 XMAX=X ( 1 ) 

SM I N = XMAX 
YM I NX= Y ( 1 ) 

YMAX=YMI NX 
VM IN = YM INX 

DO 200 1=2, NP 
I F ( X( I ) .GT.XMAX) XMAX = X( I ) 

I F ( Y ( I ) .LT.VMIN) VMI N= Y( I ) 

IF (Y( I ) .GT.YMAX) YMAX=Y(I) 

IF( X( I ) .GE.SMIN) GO TO 200 
SM I N=X ( I ) 

YMI NX=Y ( I ) 

200 CONTINUE 

C JC=NUMEER OF ENTRIES IN TABC 

IF(JC) 400,500,400 
400 DC) 900 1 = 1, JC 

IF(XMAX.LT.TABC( 1,1) .AND .YMAX . LT .TABC ( 1,2) . AND.VMIN. 
lGT.TABCt I ,3) . AND. SMI N. GT .TABCtI ,4) ) GO TO 700 
900 CONTINUE 

C DID NOT FIND THIS CONTOUR TO BE INTERIOR TO ANOTHER 

C CHECK IF EXTERIOR 

DO 1000 1=1, JC 

IF(XMAX.GT.TABC( 1,1) .AND .YMAX. GT .TABC ( I ,2) .AND.VMIN. 
1LT.TABC( I ,3) .AND.SMIN.LT.TABC( I ,4) ) GO TO 800 
1000 CONTINUE 
500 IF ( JC . EQ* 20 ) RETURN 
JC=JC+1 
MC = JC 

600 TA BC ( MC , 1 ) =XMAX 
TABC (MC ,2) =YMAX 
T A BC ( M C , 3 ) =V M I N 
TABC ( MC ,4) =S MI N 
TABC ( MC , 5 ) =YMI NX 
TABC (MC,6 ) =CV 
RETURN 

C CHECK IF THIS INTERIOR ONE IS OF HIGHER LEVEL 

700 IF(CV.LE.T ABC (1,6) ) RETURN 
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2000 

C 

800 



NIC = I 

GO TO 600 

CHECK IF LEVEL OF TH 
IF(CV.LT.TABC( I ,6) ) 
GO TO 2000 
END 



IS EXTERIOR ONE I 
RETURN 



S HIGHFR 
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